! $Id: magnetic.f90 13608 2010-04-09 11:56:42Z mppiyali $
!
! This modules deals with all aspects of magnetic fields; if no
! magnetic fields are invoked, a corresponding replacement dummy
! routine is used instead which absorbs all the calls to the
! magnetically relevant subroutines listed in here.
!
!** 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 :: lmagnetic = .true.
!
! MVAR CONTRIBUTION 3
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED aa(3); a2; aij(3,3); bb(3); bbb(3); ab; ua; uxb(3); exa(3)
! PENCILS PROVIDED b2; bij(3,3); del2a(3); graddiva(3);jj(3)
! PENCILS PROVIDED j2; jb; va2; jxb(3); jxbr(3); jxbr2; ub; uxb(3); uxb2
! PENCILS PROVIDED uxj(3); beta; djuidjbi; jo
! PENCILS PROVIDED ujxb; oxu(3); oxuxb(3); jxbxb(3); jxbrxb(3)
! PENCILS PROVIDED glnrhoxb(3); oxj(3); diva
! PENCILS PROVIDED jij(3,3); sj; ss12
! PENCILS PROVIDED etava, etaj, etaj2, etajrho
! PENCILS PROVIDED cosjb,jparallel;jperp
! PENCILS PROVIDED cosub
!
!***************************************************************
module Magnetic
!
use Cdata
use Cparam
use Messages, only: fatal_error,inevitably_fatal_error,warning,svn_id,timing
use Sub, only: keep_compiler_quiet
!
implicit none
!
include 'record_types.h'
include 'magnetic.h'
!
! Slice precalculation buffers
!
real, target, dimension (nx,ny,3) :: bb_xy, jj_xy
real, target, dimension (nx,ny,3) :: bb_xy2, jj_xy2
real, target, dimension (nx,nz,3) :: bb_xz, jj_xz
real, target, dimension (ny,nz,3) :: bb_yz, jj_yz
!
real, target, dimension (nx,ny) :: b2_xy, jb_xy
real, target, dimension (nx,ny) :: b2_xy2, jb_xy2
real, target, dimension (ny,nz) :: b2_yz, jb_yz
real, target, dimension (nx,nz) :: b2_xz, jb_xz
!
real, target, dimension (nx,ny) :: beta_xy
real, target, dimension (nx,ny) :: beta_xy2
real, target, dimension (ny,nz) :: beta_yz
real, target, dimension (nx,nz) :: beta_xz
!
! xy-averaged field
!
real, dimension (mz,3) :: aamz
real, dimension (nz,3) :: bbmz,jjmz
!
! Parameters
!
integer, parameter :: nresi_max=4
!
real, dimension (ninit) :: amplaa=0.0, kx_aa=1.0, ky_aa=1.0, kz_aa=1.0
real, dimension (ninit) :: phasex_aa=0.0, phasey_aa=0.0, phasez_aa=0.0
character (len=labellen), dimension(ninit) :: initaa='nothing'
character (len=labellen), dimension(nresi_max) :: iresistivity=''
character (len=labellen) :: fring_profile='tanh'
!
! Input parameters
!
complex, dimension(3) :: coefaa=(/0.0,0.0,0.0/), coefbb=(/0.0,0.0,0.0/)
real, dimension(3) :: B_ext=(/0.0,0.0,0.0/), B1_ext, B_ext_inv, B_ext_tmp
real, dimension(3) :: axisr1=(/0,0,1/), dispr1=(/0.0,0.5,0.0/)
real, dimension(3) :: axisr2=(/1,0,0/), dispr2=(/0.0,-0.5,0.0/)
real, dimension(3) :: axisr3=(/1,0,0/), dispr3=(/0.0,-0.5,0.0/)
real, dimension(nx,3) :: uxbb
real, target :: zmode=1.0 !(temporary)
real :: fring1=0.0, Iring1=0.0, Rring1=1.0, wr1=0.3
real :: fring2=0.0, Iring2=0.0, Rring2=1.0, wr2=0.3
real :: fring3=0.0, Iring3=0.0, Rring3=1.0, wr3=0.3
real :: radius=0.1, epsilonaa=0.01, widthaa=0.5, x0aa=0.0, z0aa=0.0
real :: by_left=0.0, by_right=0.0, bz_left=0.0, bz_right=0.0
real :: bthresh=0.0, bthresh_per_brms=0.0, brms=0.0, bthresh_scl=1.0
real :: eta_shock=0.0
real :: eta_va=0., eta_j=0., eta_j2=0., eta_jrho=0., eta_min=0., etaj20=0.
real :: rhomin_jxb=0.0, va2max_jxb=0.0
real :: omega_Bz_ext=0.0
real :: mu_r=-0.5 !(still needed for backwards compatibility)
real :: mu_ext_pot=-0.5,inclaa=0.0
real :: mu012=0.5 !(=1/2mu0)
real :: rescale_aa=0.0
real :: ampl_B0=0.0, D_smag=0.17, B_ext21, B_ext11
real :: Omega_ampl=0.0
real :: rmode=1.0, rm_int=0.0, rm_ext=0.0
real :: nu_ni=0.0, nu_ni1
real :: displacement_gun=0.0
real :: initpower_aa=0.0, cutoff_aa=0.0, brms_target=1.0
real :: rescaling_fraction=1.0
real :: phase_beltrami=0.0, ampl_beltrami=0.0
real :: bmz=0, bmz_beltrami_phase=0.0
real :: taareset=0.0, daareset=0.0
real :: center1_x=0.0, center1_y=0.0, center1_z=0.0
real :: fluxtube_border_width=impossible
real :: eta_jump=0.0
integer :: nbvec,nbvecmax=nx*ny*nz/4, va2power_jxb=5
integer :: N_modes_aa=1, naareset
integer :: nrings=2
logical :: lpress_equil=.false., lpress_equil_via_ss=.false.
logical :: lpress_equil_alt=.false.
logical :: llorentzforce=.true., linduction=.true.
logical :: lresi_eta_const=.false.
logical :: lresi_etaSS=.false.
logical :: lresi_eta_shock=.false.
logical :: lresi_etava=.false.
logical :: lresi_etaj=.false.
logical :: lresi_etaj2=.false.
logical :: lresi_etajrho=.false.
logical :: lresi_anomalous=.false.
logical, target, dimension (3) :: lfrozen_bb_bot=(/.false.,.false.,.false./)
logical, target, dimension (3) :: lfrozen_bb_top=(/.false.,.false.,.false./)
logical :: lohmic_heat=.true., lneutralion_heat=.true.
logical :: reinitialize_aa=.false.
logical :: lB_ext_pot=.false.
logical :: lforce_free_test=.false.
logical :: lgauss=.false.
logical :: lbb_as_aux=.false., ljj_as_aux=.false.
logical :: lbbt_as_aux=.false., ljjt_as_aux=.false.
logical :: lcheck_positive_va2=.false.
logical :: lreset_aa=.false.
!
namelist /magnetic_init_pars/ &
B_ext, lohmic_heat, fring1, Iring1, Rring1, wr1, axisr1, dispr1, fring2, &
Iring2, Rring2, wr2, axisr2, dispr2, fring3, Iring3, Rring3, wr3, &
axisr3, dispr3, fring_profile, nrings, radius, epsilonaa, x0aa, z0aa, &
widthaa, by_left, by_right, bz_left, bz_right, initaa, amplaa, kx_aa, &
ky_aa, kz_aa, coefaa, coefbb, phasex_aa, phasey_aa, phasez_aa, inclaa, &
lpress_equil, lpress_equil_via_ss, mu_r, mu_ext_pot, lB_ext_pot, &
lforce_free_test, ampl_B0, initpower_aa, cutoff_aa, N_modes_aa, rmode, &
zmode, rm_int, rm_ext, lgauss, lcheck_positive_va2, lbb_as_aux, &
ljj_as_aux, lbbt_as_aux, ljjt_as_aux, &
lneutralion_heat, center1_x, center1_y, center1_z, &
fluxtube_border_width, va2max_jxb, va2power_jxb, eta_jump,&
lpress_equil_alt
!
! Run parameters
!
real :: eta=0.0, eta1=0.0, eta_hyper2=0.0, eta_hyper3=0.0, eta_anom=0.0
real :: eta_int=0.0, eta_ext=0.0, wresistivity=0.01, eta_xy_max=1.0
real :: height_eta=0.0, eta_out=0.0
real :: z_surface=0.0
real :: tau_aa_exterior=0.0
real :: sigma_ratio=1.0, eta_width=0.0, eta_z0=1.0, eta_z1=1.0
real :: alphaSSm=0.0
real :: alpha_rmax=0.0, alpha_width=0.0
real :: Omega_rmax=0.0, Omega_rwidth=0.0
real :: k1_ff=1.0, ampl_ff=1.0, swirl=1.0
real :: k1x_ff=1.0, k1y_ff=1.0, k1z_ff=1.0
real :: inertial_length=0.0, linertial_2
real :: vcrit_anom=1.0
real, dimension(mx,my) :: eta_xy
real, dimension(mx,my,3) :: geta_xy
real, dimension(mz) :: coskz,sinkz,eta_z
real, dimension(mz,3) :: geta_z
logical :: lcalc_aamean=.false.
logical :: luse_Bext_in_b2=.false.
logical :: lmean_friction=.false.
logical :: llarge_scale_velocity=.false.
logical :: lhalox=.false.
logical :: lrun_initaa=.false.
character (len=labellen) :: zdep_profile='fs'
character (len=labellen) :: eta_xy_profile='schnack89'
!
namelist /magnetic_run_pars/ &
eta, eta1, eta_anom, B_ext, omega_Bz_ext, nu_ni, &
lohmic_heat,eta_out,z_surface, &
tau_aa_exterior, kx_aa, ky_aa, kz_aa, &
lcalc_aamean, k1_ff, &
ampl_ff, swirl, radius, k1x_ff, k1y_ff, k1z_ff, lcheck_positive_va2, &
lmean_friction, bthresh, bthresh_per_brms, iresistivity, &
alphaSSm, alpha_rmax, alpha_width, eta_int, &
eta_ext, eta_shock, eta_va,eta_j, eta_j2, eta_jrho, eta_min, &
wresistivity, eta_xy_max, rhomin_jxb, va2max_jxb, &
va2power_jxb, llorentzforce, linduction, reinitialize_aa, rescale_aa, &
lB_ext_pot, displacement_gun, D_smag, brms_target, &
rescaling_fraction, sigma_ratio, zdep_profile, eta_width, &
eta_z0, eta_z1, inertial_length, &
lbb_as_aux, ljj_as_aux, &
lbbt_as_aux, ljjt_as_aux, lneutralion_heat, lreset_aa, daareset, &
luse_Bext_in_b2, llarge_scale_velocity, &
lhalox, vcrit_anom, eta_jump,&
Omega_rmax,Omega_rwidth,lrun_initaa
!
! Diagnostic variables (need to be consistent with reset list below)
!
integer :: idiag_ab_int=0 ! DIAG_DOC: $\int\Av\cdot\Bv\;dV$
integer :: idiag_jb_int=0 ! DIAG_DOC: $\int\jv\cdot\Bv\;dV$
integer :: idiag_b2tm=0 ! DIAG_DOC: $\left<\bv(t)\cdot\int_0^t\bv(t')
! DIAG_DOC: dt'\right>$
integer :: idiag_bjtm=0 ! DIAG_DOC: $\left<\bv(t)\cdot\int_0^t\jv(t')
! DIAG_DOC: dt'\right>$
integer :: idiag_jbtm=0 ! DIAG_DOC: $\left<\jv(t)\cdot\int_0^t\bv(t')
! DIAG_DOC: dt'\right>$
integer :: idiag_b2ruzm=0 ! DIAG_DOC: $\left<\Bv^2\rho u_z\right>$
integer :: idiag_b2uzm=0 ! DIAG_DOC: $\left<\Bv^2u_z\right>$
integer :: idiag_ubbzm=0 ! DIAG_DOC: $\left<(\uv\cdot\Bv)B_z\right>$
integer :: idiag_b1m=0 ! DIAG_DOC: $\left<|\Bv|\right>$
integer :: idiag_b2m=0 ! DIAG_DOC: $\left<\Bv^2\right>$
integer :: idiag_bm2=0 ! DIAG_DOC: $\max(\Bv^2)$
integer :: idiag_j2m=0 ! DIAG_DOC: $\left<\jv^2\right>$
integer :: idiag_jm2=0 ! DIAG_DOC: $\max(\jv^2)$
integer :: idiag_abm=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$
integer :: idiag_abumx=0 ! DIAG_DOC: $\left<u_x\Av\cdot\Bv\right>$
integer :: idiag_abumy=0 ! DIAG_DOC: $\left<u_y\Av\cdot\Bv\right>$
integer :: idiag_abumz=0 ! DIAG_DOC: $\left<u_z\Av\cdot\Bv\right>$
integer :: idiag_abmh=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (temp)
integer :: idiag_abmn=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (north)
integer :: idiag_abms=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (south)
integer :: idiag_abrms=0 ! DIAG_DOC: $\left<(\Av\cdot\Bv)^2\right>^{1/2}$
integer :: idiag_ajm=0 ! DIAG_DOC: $\left<\jv\cdot\Av\right>$
integer :: idiag_jbm=0 ! DIAG_DOC: $\left<\jv\cdot\Bv\right>$
integer :: idiag_jbmh=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (temp)
integer :: idiag_jbmn=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (north)
integer :: idiag_jbms=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>$ (south)
integer :: idiag_ubm=0 ! DIAG_DOC: $\left<\uv\cdot\Bv\right>$
integer :: idiag_uxbxm=0 ! DIAG_DOC: $\left<u_xB_x\right>$
integer :: idiag_uybym=0 ! DIAG_DOC: $\left<u_yB_y\right>$
integer :: idiag_uzbzm=0 ! DIAG_DOC: $\left<u_zB_z\right>$
integer :: idiag_cosubm=0 ! DIAG_DOC: $\left<\Uv\cdot\Bv/(|\Uv|\,|\Bv|)\right>$
integer :: idiag_uam=0 ! DIAG_DOC: $\left<\uv\cdot\Av\right>$
integer :: idiag_ujm=0 ! DIAG_DOC: $\left<\uv\cdot\Jv\right>$
integer :: idiag_fbm=0 ! DIAG_DOC: $\left<\fv\cdot\Bv\right>$
integer :: idiag_fxbxm=0 ! DIAG_DOC: $\left<f_x B_x\right>$
integer :: idiag_epsM=0 ! DIAG_DOC: $\left<2\eta\mu_0\jv^2\right>$
integer :: idiag_epsAD=0 ! DIAG_DOC: $\left<\rho^{-1} t_{\rm AD} (\vec{J}\times\vec{B})^2\right>$ (heating by ion-neutrals friction)
integer :: idiag_bxpt=0 ! DIAG_DOC: $B_x(x_0,y_0,z_0,t)$
integer :: idiag_bypt=0 ! DIAG_DOC: $B_y(x_0,y_0,z_0,t)$
integer :: idiag_bzpt=0 ! DIAG_DOC: $B_z(x_0,y_0,z_0,t)$
integer :: idiag_Expt=0 ! DIAG_DOC: ${\cal E}_x(x_0,y_0,z_0,t)$
integer :: idiag_Eypt=0 ! DIAG_DOC: ${\cal E}_y(x_0,y_0,z_0,t)$
integer :: idiag_Ezpt=0 ! DIAG_DOC: ${\cal E}_z(x_0,y_0,z_0,t)$
integer :: idiag_epsM_LES=0 ! DIAG_DOC:
integer :: idiag_aybym2=0 ! DIAG_DOC:
integer :: idiag_exaym2=0 ! DIAG_DOC:
integer :: idiag_exabot=0 ! DIAG_DOC: $\int\Ev\times\Av\,dS|_{\rm bot}$
integer :: idiag_exatop=0 ! DIAG_DOC: $\int\Ev\times\Av\,dS|_{\rm top}$
integer :: idiag_exjm2=0 ! DIAG_DOC:
integer :: idiag_emag=0 ! DIAG_DOC: $\int_V{1\over2\mu_0}\Bv^2\, dV$
integer :: idiag_brms=0 ! DIAG_DOC: $\left<\Bv^2\right>^{1/2}$
integer :: idiag_bmax=0 ! DIAG_DOC: $\max(|\Bv|)$
integer :: idiag_bxmin=0 ! DIAG_DOC: $\min(|B_x|)$
integer :: idiag_bymin=0 ! DIAG_DOC: $\min(|B_y|)$
integer :: idiag_bzmin=0 ! DIAG_DOC: $\min(|B_z|)$
integer :: idiag_bxmax=0 ! DIAG_DOC: $\max(|B_x|)$
integer :: idiag_bymax=0 ! DIAG_DOC: $\max(|B_y|)$
integer :: idiag_bzmax=0 ! DIAG_DOC: $\max(|B_z|)$
integer :: idiag_jrms=0 ! DIAG_DOC: $\left<\jv^2\right>^{1/2}$
integer :: idiag_jmax=0 ! DIAG_DOC: $\max(|\jv|)$
integer :: idiag_vArms=0 ! DIAG_DOC: $\left<\Bv^2/\varrho\right>^{1/2}$
integer :: idiag_vAmax=0 ! DIAG_DOC: $\max(\Bv^2/\varrho)^{1/2}$
integer :: idiag_dtb=0 ! DIAG_DOC: $\delta t / [c_{\delta t}\,\delta x
! DIAG_DOC: /v_{\rm A,max}]$
! DIAG_DOC: \quad(time step relative to
! DIAG_DOC: Alfv{\'e}n time step;
! DIAG_DOC: see \S~\ref{time-step})
integer :: idiag_dteta=0 ! DIAG_DOC: $\delta t/[c_{\delta t,{\rm v}}\,
! DIAG_DOC: \delta x^2/\eta_{\rm max}]$
! DIAG_DOC: \quad(time step relative to
! DIAG_DOC: resistive time step;
! DIAG_DOC: see \S~\ref{time-step})
integer :: idiag_axm=0 ! DIAG_DOC:
integer :: idiag_aym=0 ! DIAG_DOC:
integer :: idiag_azm=0 ! DIAG_DOC:
integer :: idiag_arms=0 ! DIAG_DOC:
integer :: idiag_amax=0 ! DIAG_DOC:
integer :: idiag_Qsm=0 ! DIAG_DOC: $\left<Q_p(\overline{B})\right>$
integer :: idiag_Qpm=0 ! DIAG_DOC: $\left<Q_p(\overline{B})\right>$
integer :: idiag_qem=0 ! DIAG_DOC: $\left<q_e(\overline{B})\right>$
integer :: idiag_beta1m=0 ! DIAG_DOC: $\left<\Bv^2/(2\mu_0 p)\right>$
! DIAG_DOC: \quad(mean inverse plasma beta)
integer :: idiag_beta1max=0 ! DIAG_DOC: $\max[\Bv^2/(2\mu_0 p)]$
! DIAG_DOC: \quad(maximum inverse plasma beta)
integer :: idiag_bxm=0 ! DIAG_DOC: $\left<\left<B\right>_{yz}^2\right>_x^{1/2}$
integer :: idiag_bym=0 ! DIAG_DOC: $\left<\left<B\right>_{zx}^2\right>_y^{1/2}$
integer :: idiag_bzm=0 ! DIAG_DOC: $\left<\left<B\right>_{xy}^2\right>_z^{1/2}$
integer :: idiag_bx2m=0 ! DIAG_DOC: $\left<B_x^2\right>$
integer :: idiag_by2m=0 ! DIAG_DOC: $\left<B_y^2\right>$
integer :: idiag_bz2m=0 ! DIAG_DOC: $\left<B_z^2\right>$
integer :: idiag_bx2mx=0 ! DIAG_DOC: $\left<B_x^2\right>_{yz}$
integer :: idiag_by2mx=0 ! DIAG_DOC: $\left<B_y^2\right>_{yz}$
integer :: idiag_bz2mx=0 ! DIAG_DOC: $\left<B_z^2\right>_{yz}$
integer :: idiag_bxbym=0 ! DIAG_DOC: $\left<B_x B_y\right>$
integer :: idiag_bxbymx=0 ! DIAG_DOC: $\left<B_x B_y\right>_{yz}$
integer :: idiag_bxbzm=0 ! DIAG_DOC:
integer :: idiag_bybzm=0 ! DIAG_DOC:
integer :: idiag_djuidjbim=0 ! DIAG_DOC:
integer :: idiag_bxbymy=0 ! DIAG_DOC:
integer :: idiag_bxbzmy=0 ! DIAG_DOC:
integer :: idiag_bybzmy=0 ! DIAG_DOC:
integer :: idiag_bxbymz=0 ! DIAG_DOC:
integer :: idiag_bxbzmz=0 ! DIAG_DOC:
integer :: idiag_bybzmz=0 ! DIAG_DOC:
integer :: idiag_b2mz=0 ! DIAG_DOC:
integer :: idiag_axmz=0 ! DIAG_DOC: $\left<{\cal A}_x\right>_{xy}$
integer :: idiag_aymz=0 ! DIAG_DOC: $\left<{\cal A}_y\right>_{xy}$
integer :: idiag_azmz=0 ! DIAG_DOC: $\left<{\cal A}_z\right>_{xy}$
integer :: idiag_abuxmz=0 ! DIAG_DOC: $\left<(\Av \cdot \Bv) u_x \right>_{xy}$
integer :: idiag_abuymz=0 ! DIAG_DOC: $\left<(\Av \cdot \Bv) u_y \right>_{xy}$
integer :: idiag_abuzmz=0 ! DIAG_DOC: $\left<(\Av \cdot \Bv) u_z \right>_{xy}$
integer :: idiag_uabxmz=0 ! DIAG_DOC: $\left<(\uv \cdot \Av) B_x \right>_{xy}$
integer :: idiag_uabymz=0 ! DIAG_DOC: $\left<(\uv \cdot \Av) B_y \right>_{xy}$
integer :: idiag_uabzmz=0 ! DIAG_DOC: $\left<(\uv \cdot \Av) B_z \right>_{xy}$
integer :: idiag_bxmz=0 ! DIAG_DOC: $\left<{\cal B}_x\right>_{xy}$
integer :: idiag_bymz=0 ! DIAG_DOC: $\left<{\cal B}_y\right>_{xy}$
integer :: idiag_bzmz=0 ! DIAG_DOC: $\left<{\cal B}_z\right>_{xy}$
integer :: idiag_jxmz=0 ! DIAG_DOC: $\left<{\cal J}_x\right>_{xy}$
integer :: idiag_jymz=0 ! DIAG_DOC: $\left<{\cal J}_y\right>_{xy}$
integer :: idiag_jzmz=0 ! DIAG_DOC: $\left<{\cal J}_z\right>_{xy}$
integer :: idiag_Exmz=0 ! DIAG_DOC: $\left<{\cal E}_x\right>_{xy}$
integer :: idiag_Eymz=0 ! DIAG_DOC: $\left<{\cal E}_y\right>_{xy}$
integer :: idiag_Ezmz=0 ! DIAG_DOC: $\left<{\cal E}_z\right>_{xy}$
integer :: idiag_bmx=0 ! DIAG_DOC: $\left<\left<\Bv\right>_{yz}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $yz$-averaged
! DIAG_DOC: mean field)
integer :: idiag_bmy=0 ! DIAG_DOC: $\left<\left<\Bv\right>_{xz}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $xz$-averaged
! DIAG_DOC: mean field)
integer :: idiag_bmz=0 ! DIAG_DOC: $\left<\left<\Bv\right>_{xy}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $xy$-averaged
! DIAG_DOC: mean field)
integer :: idiag_bmzS2=0 ! DIAG_DOC: $\left<\left<\Bv_S\right>_{xy}^2\right>$
integer :: idiag_bmzA2=0 ! DIAG_DOC: $\left<\left<\Bv_A\right>_{xy}^2\right>$
integer :: idiag_jmx=0 ! DIAG_DOC: $\left<\left<\Jv\right>_{yz}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $yz$-averaged
! DIAG_DOC: mean current density)
integer :: idiag_jmy=0 ! DIAG_DOC: $\left<\left<\Jv\right>_{xz}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $xz$-averaged
! DIAG_DOC: mean current density)
integer :: idiag_jmz=0 ! DIAG_DOC: $\left<\left<\Jv\right>_{xy}^2
! DIAG_DOC: \right>^{1/2}$
! DIAG_DOC: \quad(energy of $xy$-averaged
! DIAG_DOC: mean current density)
integer :: idiag_bmzph=0 ! DIAG_DOC: Phase of a Beltrami field
integer :: idiag_bmzphe=0 ! DIAG_DOC: Error of phase of a Beltrami field
integer :: idiag_bsinphz=0 ! DIAG_DOC: sine of phase of a Beltrami field
integer :: idiag_bcosphz=0 ! DIAG_DOC: cosine of phase of a Beltrami field
integer :: idiag_emxamz3=0 ! DIAG_DOC: $\left<\left<\Ev\right>_{xy}\times\left<\Av\right>_{xy}
! DIAG_DOC: \right>$ \quad($xy$-averaged
! DIAG_DOC: mean field helicity flux)
integer :: idiag_embmz=0 ! DIAG_DOC: $\left<\left<\Ev\right>_{xy}\cdot\left<\Bv\right>_{xy}
! DIAG_DOC: \right>$ \quad($xy$-averaged
! DIAG_DOC: mean field helicity production )
integer :: idiag_ambmz=0 ! DIAG_DOC: $\left<\left<\Av\right>_{xy}\cdot\left<\Bv\right>_{xy}\right>$
! DIAG_DOC: \quad (magnetic helicity of $xy$-averaged mean field)
integer :: idiag_ambmzh=0 ! DIAG_DOC: $\left<\left<\Av\right>_{xy}\cdot\left<\Bv\right>_{xy}\right>$
! DIAG_DOC: \quad (magnetic helicity of $xy$-averaged mean field, temp)
integer :: idiag_ambmzn=0 ! DIAG_DOC: $\left<\left<\Av\right>_{xy}\cdot\left<\Bv\right>_{xy}\right>$
! DIAG_DOC: \quad (magnetic helicity of $xy$-averaged mean field, north)
integer :: idiag_ambmzs=0 ! DIAG_DOC: $\left<\left<\Av\right>_{xy}\cdot\left<\Bv\right>_{xy}\right>$
! DIAG_DOC: \quad (magnetic helicity of $xy$-averaged mean field, south)
integer :: idiag_jmbmz=0 ! DIAG_DOC: $\left<\left<\Jv\right>_{xy}\cdot\left<\Bv\right>_{xy}
! DIAG_DOC: \right>$ \quad(current helicity
! DIAG_DOC: of $xy$-averaged mean field)
integer :: idiag_kx_aa=0 ! DIAG_DOC: $k_x$
integer :: idiag_kmz=0 ! DIAG_DOC: $\left<\left<\Jv\cdot\Bv\right>_{xy}\right>/
! DIAG_DOC: \left<\left<\Jv\cdot\Bv\right>_{xy}\right>$
integer :: idiag_bx2my=0 ! DIAG_DOC: $\left< B_x^2 \right>_{xz}$
integer :: idiag_by2my=0 ! DIAG_DOC: $\left< B_y^2 \right>_{xz}$
integer :: idiag_bz2my=0 ! DIAG_DOC: $\left< B_z^2 \right>_{xz}$
integer :: idiag_bx2mz=0 ! DIAG_DOC: $\left< B_x^2 \right>_{xy}$
integer :: idiag_by2mz=0 ! DIAG_DOC: $\left< B_y^2 \right>_{xy}$
integer :: idiag_bz2mz=0 ! DIAG_DOC: $\left< B_z^2 \right>_{xy}$
integer :: idiag_bx2rmz=0 ! DIAG_DOC: $\left< B_x^2/\varrho \right>_{xy}$
integer :: idiag_by2rmz=0 ! DIAG_DOC: $\left< B_y^2/\varrho \right>_{xy}$
integer :: idiag_bz2rmz=0 ! DIAG_DOC: $\left< B_z^2/\varrho \right>_{xy}$
integer :: idiag_jbmz=0 ! DIAG_DOC: $\left<\Jv\cdot\Bv\right>|_{xy}$
integer :: idiag_abmz=0 ! DIAG_DOC: $\left<\Av\cdot\Bv\right>|_{xy}$
integer :: idiag_uamz=0 ! DIAG_DOC: $\left<\uv\cdot\Av\right>|_{xy}$
integer :: idiag_examz1=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>_{xy}|_x$
integer :: idiag_examz2=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>_{xy}|_y$
integer :: idiag_examz3=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>_{xy}|_z$
integer :: idiag_bxmxy=0 ! DIAG_DOC: $\left< B_x \right>_{xy}$
integer :: idiag_bymxy=0 ! DIAG_DOC: $\left< B_y \right>_{xy}$
integer :: idiag_bzmxy=0 ! DIAG_DOC: $\left< B_z \right>_{xy}$
integer :: idiag_jxmxy=0 ! DIAG_DOC: $\left< J_x \right>_{xy}$
integer :: idiag_jymxy=0 ! DIAG_DOC: $\left< J_y \right>_{xy}$
integer :: idiag_jzmxy=0 ! DIAG_DOC: $\left< J_z \right>_{xy}$
integer :: idiag_axmxz=0 ! DIAG_DOC: $\left< A_x \right>_{xz}$
integer :: idiag_aymxz=0 ! DIAG_DOC: $\left< A_y \right>_{xz}$
integer :: idiag_azmxz=0 ! DIAG_DOC: $\left< A_z \right>_{xz}$
integer :: idiag_bxmxz=0 ! DIAG_DOC: $\left< B_x \right>_{xz}$
integer :: idiag_bymxz=0 ! DIAG_DOC: $\left< B_y \right>_{xz}$
integer :: idiag_bzmxz=0 ! DIAG_DOC: $\left< B_z \right>_{xz}$
integer :: idiag_bx2mxz=0 ! DIAG_DOC: $\left< B_x^2 \right>_{xz}$
integer :: idiag_by2mxz=0 ! DIAG_DOC: $\left< B_y^2 \right>_{xz}$
integer :: idiag_bz2mxz=0 ! DIAG_DOC: $\left< B_z^2 \right>_{xz}$
integer :: idiag_bx2mxy=0 ! DIAG_DOC: $\left< B_x^2 \right>_{z}$
integer :: idiag_by2mxy=0 ! DIAG_DOC: $\left< B_y^2 \right>_{z}$
integer :: idiag_bz2mxy=0 ! DIAG_DOC: $\left< B_z^2 \right>_{z}$
integer :: idiag_bxbymxy=0 ! DIAG_DOC: $\left< B_x B_y \right>_{z}$
integer :: idiag_bxbzmxy=0 ! DIAG_DOC: $\left< B_x B_z \right>_{z}$
integer :: idiag_bybzmxy=0 ! DIAG_DOC: $\left< B_y B_z \right>_{z}$
integer :: idiag_bxbymxz=0 ! DIAG_DOC: $\left< B_x B_y \right>_{y}$
integer :: idiag_bxbzmxz=0 ! DIAG_DOC: $\left< B_x B_z \right>_{y}$
integer :: idiag_bybzmxz=0 ! DIAG_DOC: $\left< B_y B_z \right>_{y}$
integer :: idiag_uxbm=0 ! DIAG_DOC: $\left<\uv\times\Bv\right>\cdot\Bv_0/B_0^2$
integer :: idiag_jxbm=0 ! DIAG_DOC: $\left<\jv\times\Bv\right>\cdot\Bv_0/B_0^2$
integer :: idiag_oxuxbm=0 ! DIAG_DOC:
integer :: idiag_jxbxbm=0 ! DIAG_DOC:
integer :: idiag_gpxbm=0 ! DIAG_DOC:
integer :: idiag_uxDxuxbm=0 ! DIAG_DOC:
integer :: idiag_b3b21m=0 ! DIAG_DOC: $\left<B_3 B_{2,1} \right>$
integer :: idiag_b1b32m=0 ! DIAG_DOC: $\left<B_1 B_{3,2} \right>$
integer :: idiag_b2b13m=0 ! DIAG_DOC: $\left<B_2 B_{1,3} \right>$
integer :: idiag_udotxbm=0 ! DIAG_DOC:
integer :: idiag_uxbdotm=0 ! DIAG_DOC:
integer :: idiag_uxbmx=0 ! DIAG_DOC: $\left<(\uv\times\Bv)_x\right>$
integer :: idiag_uxbmy=0 ! DIAG_DOC: $\left<(\uv\times\Bv)_y\right>$
integer :: idiag_uxbmz=0 ! DIAG_DOC: $\left<(\uv\times\Bv)_z\right>$
integer :: idiag_jxbmx=0 ! DIAG_DOC: $\left<(\jv\times\Bv)_x\right>$
integer :: idiag_jxbmy=0 ! DIAG_DOC: $\left<(\jv\times\Bv)_y\right>$
integer :: idiag_jxbmz=0 ! DIAG_DOC: $\left<(\jv\times\Bv)_z\right>$
integer :: idiag_uxbcmx=0 ! DIAG_DOC:
integer :: idiag_uxbcmy=0 ! DIAG_DOC:
integer :: idiag_uxbsmx=0 ! DIAG_DOC:
integer :: idiag_uxbsmy=0 ! DIAG_DOC:
integer :: idiag_examx=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>|_x$
integer :: idiag_examy=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>|_y$
integer :: idiag_examz=0 ! DIAG_DOC: $\left<\Ev\times\Av\right>|_z$
integer :: idiag_exjmx=0 ! DIAG_DOC: $\left<\Ev\times\Jv\right>|_x$
integer :: idiag_exjmy=0 ! DIAG_DOC: $\left<\Ev\times\Jv\right>|_y$
integer :: idiag_exjmz=0 ! DIAG_DOC: $\left<\Ev\times\Jv\right>|_z$
integer :: idiag_dexbmx=0 ! DIAG_DOC: $\left<\nabla\times\Ev\times\Bv\right>|_x$
integer :: idiag_dexbmy=0 ! DIAG_DOC: $\left<\nabla\times\Ev\times\Bv\right>|_y$
integer :: idiag_dexbmz=0 ! DIAG_DOC: $\left<\nabla\times\Ev\times\Bv\right>|_z$
integer :: idiag_phibmx=0 ! DIAG_DOC: $\left<\phi\Bv\right>|_x$
integer :: idiag_phibmy=0 ! DIAG_DOC: $\left<\phi\Bv\right>|_y$
integer :: idiag_phibmz=0 ! DIAG_DOC: $\left<\phi\Bv\right>|_z$
integer :: idiag_uxjm=0 ! DIAG_DOC:
integer :: idiag_b2divum=0 ! DIAG_DOC: $\left<\Bv^2\nabla\cdot\uv\right>$
integer :: idiag_ujxbm=0 ! DIAG_DOC: $\left<\uv\cdot(\Jv\times\Bv\right>$
integer :: idiag_jxbrxm=0 ! DIAG_DOC:
integer :: idiag_jxbrym=0 ! DIAG_DOC:
integer :: idiag_jxbrzm=0 ! DIAG_DOC:
integer :: idiag_jxbr2m=0 ! DIAG_DOC: $\left<(\Jv\times\Bv/\rho)^2\right>$
integer :: idiag_jxbrxmx=0 ! DIAG_DOC:
integer :: idiag_jxbrymx=0 ! DIAG_DOC:
integer :: idiag_jxbrzmx=0 ! DIAG_DOC:
integer :: idiag_jxbrxmy=0 ! DIAG_DOC:
integer :: idiag_jxbrymy=0 ! DIAG_DOC:
integer :: idiag_jxbrzmy=0 ! DIAG_DOC:
integer :: idiag_jxbrxmz=0 ! DIAG_DOC:
integer :: idiag_jxbrymz=0 ! DIAG_DOC:
integer :: idiag_jxbrzmz=0 ! DIAG_DOC:
integer :: idiag_uxBrms=0 ! DIAG_DOC:
integer :: idiag_Bresrms=0 ! DIAG_DOC:
integer :: idiag_Rmrms=0 ! DIAG_DOC:
integer :: idiag_jfm=0 ! DIAG_DOC:
integer :: idiag_vA2m=0 ! DIAG_DOC:
integer :: idiag_bxmx=0 ! DIAG_DOC:
integer :: idiag_bymx=0 ! DIAG_DOC:
integer :: idiag_bzmx=0 ! DIAG_DOC:
integer :: idiag_bxmy=0 ! DIAG_DOC:
integer :: idiag_bymy=0 ! DIAG_DOC:
integer :: idiag_bzmy=0 ! DIAG_DOC:
integer :: idiag_mflux_x=0 ! DIAG_DOC:
integer :: idiag_mflux_y=0 ! DIAG_DOC:
integer :: idiag_mflux_z=0 ! DIAG_DOC:
integer :: idiag_bmxy_rms=0 ! DIAG_DOC: $\sqrt{[\left<b_x\right>_z(x,y)]^2 +
! DIAG_DOC: [\left<b_y\right>_z(x,y)]^2 +
! DIAG_DOC: [\left<b_z\right>_z(x,y)]^2} $
integer :: idiag_etasmagm=0 ! DIAG_DOC: Mean of Smagorinsky resistivity
integer :: idiag_etasmagmin=0 ! DIAG_DOC: Min of Smagorinsky resistivity
integer :: idiag_etasmagmax=0 ! DIAG_DOC: Max of Smagorinsky resistivity
integer :: idiag_etavamax=0 ! DIAG_DOC: Max of artificial resistivity
! DIAG_DOC: $\eta\sim v_A$
integer :: idiag_etajmax=0 ! DIAG_DOC: Max of artificial resistivity
! DIAG_DOC: $\eta\sim J / \sqrt{\rho}$
integer :: idiag_etaj2max=0 ! DIAG_DOC: Max of artificial resistivity
! DIAG_DOC: $\eta\sim J^2 / \rho$
integer :: idiag_etajrhomax=0 ! DIAG_DOC: Max of artificial resistivity
! DIAG_DOC: $\eta\sim J / \rho$
integer :: idiag_cosjbm=0 ! DIAG_DOC: $\left<\Jv\cdot\Bv/(|\Jv|\,|\Bv|)\right>$
integer :: idiag_jparallelm=0 ! DIAG_DOC: Mean value of the component
! DIAG_DOC: of J parallel to B
integer :: idiag_jperpm=0 ! DIAG_DOC: Mean value of the component
! DIAG_DOC: of J perpendicular to B
integer :: idiag_brmsn=0,idiag_brmss=0,idiag_brmsh=0
integer :: idiag_Exmxz=0 ! DIAG_DOC: $\left<{\cal E}_x\right>_{y}$
integer :: idiag_Eymxz=0 ! DIAG_DOC: $\left<{\cal E}_y\right>_{y}$
integer :: idiag_Ezmxz=0 ! DIAG_DOC: $\left<{\cal E}_z\right>_{y}$
integer :: idiag_etatotalmx=0 ! DIAG_DOC: $\left<\eta\right>_{yz}$
integer :: idiag_etatotalmz=0 ! DIAG_DOC: $\left<\eta\right>_{xy}$
!
contains
!***********************************************************************
subroutine register_magnetic()
!
! Initialise variables which should know that we solve for the vector
! potential: iaa, etc; increase nvar accordingly
!
! 1-may-02/wolf: coded
!
use FArrayManager, only: farray_register_pde
!
call farray_register_pde('aa',iaa,vector=3)
iax = iaa; iay = iaa+1; iaz = iaa+2
!
! Identify version number.
!
if (lroot) call svn_id( &
"$Id: magnetic.f90 13608 2010-04-09 11:56:42Z mppiyali $")
!
! Writing files for use with IDL
!
if (lroot) then
if (maux == 0) then
if (nvar < mvar) write(4,*) ',aa $'
if (nvar == mvar) write(4,*) ',aa'
else
write(4,*) ',aa $'
endif
write(15,*) 'aa = fltarr(mx,my,mz,3)*one'
endif
!
endsubroutine register_magnetic
!***********************************************************************
subroutine initialize_magnetic(f,lstarting)
!
! Perform any post-parameter-read initialization
!
! 24-nov-02/tony: dummy routine - nothing to do at present
! 20-may-03/axel: reinitialize_aa added
!
use FArrayManager
use SharedVariables
use EquationOfState, only: cs0
!
real, dimension (mx,my,mz,mfarray) :: f
logical :: lstarting
integer :: i
!
! Precalculate 1/mu (moved here from register.f90)
!
mu01=1./mu0
mu012=.5*mu01
!
! Precalculate eta if 1/eta (==eta1) is given instead
!
if (eta1/=0.0) then
eta=1./eta1
endif
!
! Precalculate 1/inertial_length^2
!
if (inertial_length/=0.0) then
linertial_2 = inertial_length**(-2)
else
linertial_2 = 0.0
! make sure not to use this value by checking that
! (inertial_length /= 0.)...
endif
!
! Precalculate 1/nu_ni
!
if (nu_ni/=0.0) then
nu_ni1=1./nu_ni
else
nu_ni1=0.0
endif
!
! calculate B_ext21 = 1/B_ext**2 and the unit vector B1_ext = B_ext/|B_ext|
! Also calculate B_ext_inv = B_ext/|B_ext|^2
!
B_ext21=B_ext(1)**2+B_ext(2)**2+B_ext(3)**2
if (B_ext21/=0.0) then
B_ext21=1/B_ext21
else
B_ext21=1.0
endif
B_ext11=sqrt(B_ext21)
B1_ext=B_ext*B_ext11
B_ext_inv=B_ext*B_ext21
!
! Initialize resistivity.
!
if (iresistivity(1)=='') iresistivity(1)='eta-const' ! default
lresi_eta_const=.false.
lresi_eta_shock=.false.
!
do i=1,nresi_max
select case (iresistivity(i))
case ('eta-const')
if (lroot) print*, 'resistivity: constant eta'
lresi_eta_const=.true.
case ('shock','eta-shock')
if (lroot) print*, 'resistivity: shock'
lresi_eta_shock=.true.
if (.not. lshock) &
call fatal_error('initialize_magnetic', &
'shock resistivity, but module setting SHOCK=noshock')
case ('none')
! do nothing
case ('')
! do nothing
case default
if (lroot) print*, 'No such such value for iresistivity(',i,'): ', &
trim(iresistivity(i))
call fatal_error('initialize_magnetic','')
endselect
enddo
!
! If we're timestepping, die or warn if the the resistivity coefficient that
! corresponds to the chosen resistivity type is not set.
!
if (lrun) then
if (lresi_eta_const.and.(eta==0.0)) &
call warning('initialize_magnetic', &
'Resistivity coefficient eta is zero!')
if (lresi_eta_shock.and.eta_shock==0.0) &
call fatal_error('initialize_magnetic', &
'Resistivity coefficient eta_shock is zero!')
endif
!
! Register an extra aux slot for bb if requested (so bb and jj are written
! to snapshots and can be easily analyzed later). For this to work you
! must reserve enough auxiliary workspace by setting, for example,
! ! MAUX CONTRIBUTION 6
! in the beginning of your src/cparam.local file, *before* setting
! ncpus, nprocy, etc.
!
! After a reload, we need to rewrite index.pro, but the auxiliary
! arrays are already allocated and must not be allocated again.
!
if (lbb_as_aux) then
if (ibb==0) then
call farray_register_auxiliary('bb',ibb,vector=3)
ibx=ibb
iby=ibb+1
ibz=ibb+2
endif
if (ibb/=0.and.lroot) then
print*, 'initialize_magnetic: ibb = ', ibb
open(3,file=trim(datadir)//'/index.pro', POSITION='append')
write(3,*) 'ibb=',ibb
write(3,*) 'ibx=',ibx
write(3,*) 'iby=',iby
write(3,*) 'ibz=',ibz
close(3)
endif
endif
!
! do the same for jj (current density)
!
if (ljj_as_aux) then
if (ijj==0) then
call farray_register_auxiliary('jj',ijj,vector=3)
ijx=ijj
ijy=ijj+1
ijz=ijj+2
endif
if (ijj/=0.and.lroot) then
print*, 'initialize_magnetic: ijj = ', ijj
open(3,file=trim(datadir)//'/index.pro', POSITION='append')
write(3,*) 'ijj=',ijj
write(3,*) 'ijx=',ijx
write(3,*) 'ijy=',ijy
write(3,*) 'ijz=',ijz
close(3)
endif
endif
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(lstarting)
!
endsubroutine initialize_magnetic
!***********************************************************************
subroutine init_aa(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 EquationOfState
use FArrayManager
use Gravity, only: gravz, z1, z2
use Initcond
use Boundcond
use InitialCondition, only: initial_condition_aa
use Mpicomm
use SharedVariables
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
!
real, dimension (mz) :: tmp
integer :: j
!
do j=1,ninit
!
select case (initaa(j))
case ('nothing'); if (lroot .and. j==1) print*,'init_aa: nothing'
case ('zero', '0'); f(:,:,:,iax:iaz) = 0.0
case ('rescale'); f(:,:,:,iax:iaz)=amplaa(j)*f(:,:,:,iax:iaz)
case ('bsiny'); call acosy(amplaa(j),f,iaa,ky_aa(j))
case ('mode'); call modev(amplaa(j),coefaa,f,iaa,kx_aa(j),ky_aa(j),kz_aa(j))
case ('modeb'); call modeb(amplaa(j),coefbb,f,iaa,kx_aa(j),ky_aa(j),kz_aa(j))
case ('sph_constb'); call sph_constb(amplaa(j),f,iaa)
case ('random-isotropic-KS')
call random_isotropic_KS(initpower_aa,f,iax,N_modes_aa)
case ('gaussian-noise'); call gaunoise(amplaa(j),f,iax,iaz)
case ('gaussian-noise-zprof')
tmp=amplaa(1)*0.5*(tanh((z-z1)/0.05)-tanh((z-z2)/0.05))
call gaunoise(tmp,f,iax,iaz)
!
! Beltrami fields, put k=-k to make sure B=curl(A) has the right phase
!
case ('Beltrami-x')
call beltrami(amplaa(j),f,iaa,KX=-kx_aa(j),phase=phasex_aa(j))
case ('Beltrami-xy')
call beltrami(amplaa(j),f,iaa,KX=kx_aa(j),phase=phasex_aa(j))
call beltrami(amplaa(j),f,iaa,KY=-kx_aa(j),phase=phasex_aa(j))
case ('Beltrami-y'); call beltrami(amplaa(j),f,iaa,KY=-ky_aa(j),phase=phasey_aa(j))
case ('Beltrami-z'); call beltrami(amplaa(j),f,iaa,KZ=-kz_aa(j),phase=phasez_aa(j))
!
case ('propto-ux'); call wave_uu(amplaa(j),f,iaa,kx=kx_aa(j))
case ('propto-uy'); call wave_uu(amplaa(j),f,iaa,ky=ky_aa(j))
case ('propto-uz'); call wave_uu(amplaa(j),f,iaa,kz=kz_aa(j))
case ('diffrot'); call diffrot(amplaa(j),f,iay)
case ('hor-tube'); call htube(amplaa(j),f,iax,iaz,radius,epsilonaa, &
center1_x,center1_z)
case ('hor-tube-x'); call htube_x(amplaa(j),f,iax,iaz,radius,epsilonaa, &
center1_y,center1_z)
case ('hor-tube_erf'); call htube_erf(amplaa(j),f,iax,iaz,radius,epsilonaa, &
center1_x,center1_z,fluxtube_border_width)
case ('hor-fluxlayer'); call hfluxlayer(amplaa(j),f,iaa,z0aa,widthaa)
case ('hor-fluxlayer-y'); call hfluxlayer_y(amplaa(j),f,iaa,z0aa,widthaa)
case ('ver-fluxlayer'); call vfluxlayer(amplaa(j),f,iaa,x0aa,widthaa)
case ('mag-support'); call magsupport(amplaa(j),f,gravz,cs0,rho0)
case ('arcade-x'); call arcade_x(amplaa(j),f,iaa,kx_aa(j),kz_aa(j))
case ('halfcos-Bx'); call halfcos_x(amplaa(j),f,iaa)
case ('uniform-Bx'); call uniform_x(amplaa(j),f,iaa)
case ('uniform-By'); call uniform_y(amplaa(j),f,iaa)
case ('uniform-Bz'); call uniform_z(amplaa(j),f,iaa)
case ('uniform-Bphi'); call uniform_phi(amplaa(j),f,iaa)
case ('Bz(x)', '3'); call vfield(amplaa(j),f,iaa)
case ('vfield2'); call vfield2(amplaa(j),f,iaa)
case ('bipolar'); call bipolar(amplaa(j),f,iaa,kx_aa(j),ky_aa(j),kz_aa(j))
case ('bipolar_restzero'); call bipolar_restzero(amplaa(j),f,iaa,kx_aa(j),ky_aa(j))
case ('vecpatternxy'); call vecpatternxy(amplaa(j),f,iaa,kx_aa(j),ky_aa(j),kz_aa(j))
case ('xjump'); call bjump(f,iaa,by_left,by_right,bz_left,bz_right,widthaa,'x')
case ('sinxsinz'); call sinxsinz(amplaa(j),f,iaa,kx_aa(j),ky_aa(j),kz_aa(j))
case ('sinxsinz_Hz'); call sinxsinz(amplaa(j),f,iaa,kx_aa(j),ky_aa(j),kz_aa(j),KKz=kz_aa(j))
case ('sin2xsin2y'); call sin2x_sin2y_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.)
case ('cosxcosy'); call cosx_cosy_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.)
case ('sinxsiny'); call sinx_siny_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.)
case ('xsiny'); call x_siny_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.)
case ('x1siny'); call x1_siny_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.,phasey_aa(j))
case ('sinxcosz'); call sinx_siny_cosz(amplaa(j),f,iay,kx_aa(j),ky_aa(j),kz_aa(j))
case ('sinycosz'); call cosx_siny_cosz(amplaa(j),f,iax,kx_aa(j),ky_aa(j),0.)
case ('cosysinz'); call cosy_sinz(amplaa(j),f,iax,ky_aa(j),kz_aa(j))
case ('x3cosycosz'); call x3_cosy_cosz(amplaa(j),f,iax,ky_aa(j),kz_aa(j))
case ('Ax=cosysinz'); call cosy_sinz(amplaa(j),f,iax,ky_aa(j),kz_aa(j))
case ('cosxcoscosy'); call cosx_coscosy_cosz(amplaa(j),f,iaz,kx_aa(j),ky_aa(j),0.)
case ('crazy', '5'); call crazy(amplaa(j),f,iaa)
!DM
case ('strange'); call strange(amplaa(j),f,iaa)
case ('sinwave-x'); call sinwave(amplaa(j),f,iaa,kx=kx_aa(j))
case ('coswave-Ax-kx'); call coswave(amplaa(j),f,iax,kx=kx_aa(j))
case ('coswave-Ax-ky'); call coswave(amplaa(j),f,iax,ky=ky_aa(j))
case ('coswave-Ax-kz'); call coswave(amplaa(j),f,iax,kz=kz_aa(j))
case ('coswave-Ay-kx'); call coswave(amplaa(j),f,iay,kx=kx_aa(j))
case ('coswave-Ay-ky'); call coswave(amplaa(j),f,iay,ky=ky_aa(j))
case ('coswave-Ay-kz'); call coswave(amplaa(j),f,iay,kz=kz_aa(j))
case ('coswave-Az-kx'); call coswave(amplaa(j),f,iaz,kx=kx_aa(j))
case ('coswave-Az-ky'); call coswave(amplaa(j),f,iaz,ky=ky_aa(j))
case ('coswave-Az-kz'); call coswave(amplaa(j),f,iaz,kz=kz_aa(j))
case ('sinwave-Ay-kz'); call sinwave_phase(f,iay,amplaa(j),kx_aa(j),ky_aa(j),kz_aa(j),phasez_aa(j))
case ('linear-zx')
do n=n1,n2; do m=m1,m2
f(l1:l2,m,n,iay)=-0.5*amplaa(j)*z(n)**2/Lxyz(3)
enddo; enddo
case ('Ferriere-uniform-Bx'); call ferriere_uniform_x(amplaa(j),f,iaa)
case ('Ferriere-uniform-By'); call ferriere_uniform_y(amplaa(j),f,iaa)
!
case default
!
! Catch unknown values
!
call fatal_error('init_aa', &
'init_aa value "' // trim(initaa(j)) // '" not recognised')
endselect
!
! End loop over initial conditions
!
enddo
!
! Interface for user's own initial condition
!
if (linitial_condition) call initial_condition_aa(f)
!
endsubroutine init_aa
!***********************************************************************
subroutine pencil_criteria_magnetic()
!
! All pencils that the Magnetic module depends on are specified here.
!
! 19-nov-04/anders: coded
!
use Mpicomm, only: stop_it
!
lpenc_requested(i_bb)=.true.
lpenc_requested(i_uxb)=.true.
!
if (dvid/=0.0) lpenc_video(i_b2)=.true.
if (dvid/=0.0) lpenc_video(i_jb)=.true.
!
! jj pencil always needed when in Weyl gauge
!
if (ietat/=0) lpenc_requested(i_del2a)=.true.
if (lresi_eta_const.or.lresi_eta_shock) &
lpenc_requested(i_del2a)=.true.
if (lresi_eta_shock) then
lpenc_requested(i_shock)=.true.
lpenc_requested(i_gshock)=.true.
lpenc_requested(i_diva)=.true.
endif
if (lresi_etava) then
lpenc_requested(i_etava)=.true.
lpenc_requested(i_jj)=.true.
endif
if (lresi_etaj) then
lpenc_requested(i_etaj)=.true.
lpenc_requested(i_jj)=.true.
endif
if (lresi_etaj2) then
lpenc_requested(i_etaj2)=.true.
lpenc_requested(i_jj)=.true.
endif
if (lresi_etajrho) then
lpenc_requested(i_etajrho)=.true.
lpenc_requested(i_jj)=.true.
endif
if (lentropy.and.lohmic_heat) lpenc_requested(i_j2)=.true.
if (lresi_anomalous) then
lpenc_requested(i_jj)=.true.
lpenc_requested(i_rho1)=.true.
endif
if (lentropy.or.ldt) lpenc_requested(i_rho1)=.true.
if (lentropy) lpenc_requested(i_TT1)=.true.
!
! ambipolar diffusion
!
if (nu_ni/=0.0) then
lpenc_requested(i_va2)=.true.
lpenc_requested(i_jxbrxb)=.true.
lpenc_requested(i_jxbr2)=.true.
endif
if ((lhydro .and. llorentzforce) .or. nu_ni/=0.0) &
lpenc_requested(i_jxbr)=.true.
if (nu_ni/=0.0) lpenc_requested(i_va2)=.true.
!
if (idiag_jxbrxm/=0 .or. idiag_jxbrym/=0 .or. idiag_jxbrzm/=0) &
lpenc_diagnos(i_jxbr)=.true.
if (idiag_jxbr2m/=0) lpenc_diagnos(i_jxbr2)=.true.
if (idiag_jxbrxmx/=0 .or. idiag_jxbrymx/=0 .or. idiag_jxbrzmx/=0 .or. &
idiag_jxbrxmy/=0 .or. idiag_jxbrymy/=0 .or. idiag_jxbrzmy/=0 .or. &
idiag_jxbrxmz/=0 .or. idiag_jxbrymz/=0 .or. idiag_jxbrzmz/=0) &
lpenc_diagnos(i_jxbr)=.true.
!
if (idiag_b2ruzm/=0) &
lpenc_diagnos(i_rho)=.true.
!
if (idiag_axmxz/=0 .or. idiag_aymxz/=0 .or. idiag_azmxz/=0) &
lpenc_diagnos2d(i_aa)=.true.
!
if (idiag_aybym2/=0 .or. idiag_exaym2/=0 .or. &
idiag_examx/=0 .or. idiag_examy/=0 .or. idiag_examz/=0 .or. &
idiag_examz1/=0 .or. idiag_examz2/=0 .or. idiag_examz3/=0 &
) lpenc_diagnos(i_aa)=.true.
if (idiag_arms/=0 .or. idiag_amax/=0) lpenc_diagnos(i_a2)=.true.
if (idiag_ab_int/=0 .or. idiag_abm/=0 .or. idiag_abmh/=0 &
.or. idiag_abmz/=0 .or. idiag_abrms/=0 &
.or. idiag_abumx/=0 .or. idiag_abumy/=0 .or. idiag_abumz/=0 &
.or. idiag_abuxmz/=0 .or. idiag_abuymz/=0 .or. idiag_abuzmz/=0 &
) lpenc_diagnos(i_ab)=.true.
if (idiag_uam/=0 .or. idiag_uamz/=0) lpenc_diagnos(i_ua)=.true.
if (idiag_djuidjbim/=0 .or. idiag_b3b21m/=0 &
.or. idiag_dexbmx/=0 .or. idiag_dexbmy/=0 .or. idiag_dexbmz/=0 &
.or. idiag_b1b32m/=0 .or. idiag_b2b13m/=0) &
lpenc_diagnos(i_bij)=.true.
if (idiag_j2m/=0 .or. idiag_jm2/=0 .or. idiag_jrms/=0 .or. &
idiag_jmax/=0 .or. idiag_epsM/=0 .or. idiag_epsM_LES/=0 .or. &
idiag_ajm/=0 ) &
lpenc_diagnos(i_j2)=.true.
if (idiag_epsAD/=0) lpenc_diagnos(i_jxbr2)=.true.
if (idiag_jb_int/=0 .or. idiag_jbm/=0 .or. idiag_jbmz/=0) lpenc_diagnos(i_jb)=.true.
if (idiag_vArms/=0 .or. idiag_vAmax/=0 .or. idiag_vA2m/=0) lpenc_diagnos(i_va2)=.true.
if (idiag_cosubm/=0) lpenc_diagnos(i_cosub)=.true.
if (idiag_ubm/=0 .or. idiag_ubbzm/=0) lpenc_diagnos(i_ub)=.true.
if (idiag_djuidjbim/=0 .or. idiag_uxDxuxbm/=0) lpenc_diagnos(i_uij)=.true.
if (idiag_uxjm/=0) lpenc_diagnos(i_uxj)=.true.
if (idiag_vArms/=0 .or. idiag_vAmax/=0) lpenc_diagnos(i_va2)=.true.
if (idiag_uxBrms/=0 .or. idiag_Rmrms/=0) lpenc_diagnos(i_uxb2)=.true.
if (idiag_beta1m/=0 .or. idiag_beta1max/=0) lpenc_diagnos(i_beta)=.true.
if (idiag_djuidjbim/=0) lpenc_diagnos(i_djuidjbi)=.true.
if (idiag_b2divum/=0) lpenc_diagnos(i_divu)=.true.
if (idiag_b2divum/=0) lpenc_diagnos(i_b2)=.true.
if (idiag_ujxbm/=0) lpenc_diagnos(i_ujxb)=.true.
if (idiag_gpxbm/=0) lpenc_diagnos(i_glnrhoxb)=.true.
if (idiag_jxbxbm/=0) lpenc_diagnos(i_jxbxb)=.true.
if (idiag_oxuxbm/=0) lpenc_diagnos(i_oxuxb)=.true.
if (idiag_exaym2/=0 .or. idiag_exjm2/=0 &
.or. idiag_jmx/=0 .or. idiag_jmy/=0 .or. idiag_jmz/=0 &
.or. idiag_ambmz/=0 .or. idiag_jmbmz/=0 .or. idiag_kmz/=0 &
.or. idiag_examx/=0 .or. idiag_examy/=0 .or. idiag_examz/=0 &
.or. idiag_examz1/=0 .or. idiag_examz2/=0 .or. idiag_examz3/=0 &
.or. idiag_exjmx/=0 .or. idiag_exjmy/=0 .or. idiag_exjmz/=0 &
) lpenc_diagnos(i_jj)=.true.
if (idiag_examz1/=0 .or. idiag_examz2/=0 .or. idiag_examz3/=0 &
) lpenc_diagnos(i_exa)=.true.
if (idiag_phibmx/=0 .or. idiag_phibmy/=0 .or. idiag_phibmz/=0 &
) lpenc_diagnos(i_diva)=.true.
if (idiag_b2uzm/=0 .or. idiag_b2ruzm/=0 .or. &
idiag_b1m/=0 .or. idiag_b2m/=0 .or. idiag_bm2/=0 .or. &
idiag_brmsh/=0 .or. idiag_brmsn/=0 .or. idiag_brmss/=0 .or. &
idiag_brms/=0 .or. idiag_bmax/=0 .or. &
idiag_emag/=0 ) &
lpenc_diagnos(i_b2)=.true.
if (idiag_etavamax/=0) lpenc_diagnos(i_etava)=.true.
if (idiag_etajmax/=0) lpenc_diagnos(i_etaj)=.true.
if (idiag_etaj2max/=0) lpenc_diagnos(i_etaj2)=.true.
if (idiag_etajrhomax/=0) lpenc_diagnos(i_etajrho)=.true.
if (idiag_cosjbm/=0) lpenc_diagnos(i_cosjb)=.true.
if (idiag_cosubm/=0) then
lpenc_diagnos(i_cosub)=.true.
endif
if ((idiag_jparallelm/=0).or.(idiag_jperpm/=0)) then
lpenc_diagnos(i_cosjb)=.true.
lpenc_diagnos(i_jparallel)=.true.
lpenc_diagnos(i_jperp)=.true.
endif
if (idiag_abumx/=0 .or. idiag_abumy/=0 .or. idiag_abumz/=0 &
.or. idiag_abuxmz/=0 .or. idiag_abuymz/=0 .or. idiag_abuzmz/=0) &
lpenc_diagnos(i_uu)=.true.
!
! Check whether right variables are set for half-box calculations.
!
if (idiag_brmsn/=0 .or. idiag_abmn/=0 .or. idiag_ambmzn/=0 &
.or. idiag_jbmn/= 0 ) then
if ((.not.lequatory).and.(.not.lequatorz)) then
call stop_it("You have to set either of lequatory or lequatorz to true to calculate averages over half the box")
else
if (lequatory) write(*,*) 'pencil-criteria_magnetic: box divided along y dirn'
if (lequatorz) write(*,*) 'pencil-criteria_magnetic: box divided along z dirn'
endif
else
endif
!
endsubroutine pencil_criteria_magnetic
!***********************************************************************
subroutine pencil_interdep_magnetic(lpencil_in)
!
! Interdependency among pencils from the Magnetic module is specified here.
!
! 19-11-04/anders: coded
!
logical, dimension(npencils) :: lpencil_in
!
if (lpencil_in(i_jparallel).or.lpencil_in(i_jperp)) then
lpencil_in(i_cosjb)=.true.
lpencil_in(i_jxb)=.true.
endif
if (lpencil_in(i_cosjb)) then
lpencil_in(i_b2)=.true.
lpencil_in(i_j2)=.true.
lpencil_in(i_jb)=.true.
endif
if (lpencil_in(i_a2)) lpencil_in(i_aa)=.true.
if (lpencil_in(i_ab)) then
lpencil_in(i_aa)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_ua)) then
lpencil_in(i_uu)=.true.
lpencil_in(i_aa)=.true.
endif
if (lpencil_in(i_va2)) then
lpencil_in(i_b2)=.true.
lpencil_in(i_rho1)=.true.
endif
if (lpencil_in(i_etava)) lpencil_in(i_va2)=.true.
if (lpencil_in(i_etaj) .or. lpencil_in(i_etaj2) .or. lpencil_in(i_etajrho)) then
lpencil_in(i_j2)=.true.
lpencil_in(i_rho1)=.true.
endif
if (lpencil_in(i_j2)) lpencil_in(i_jj)=.true.
if (lpencil_in(i_uxj)) then
lpencil_in(i_uu)=.true.
lpencil_in(i_jj)=.true.
endif
if (lpencil_in(i_jb)) then
lpencil_in(i_bb)=.true.
lpencil_in(i_jj)=.true.
endif
if (lpencil_in(i_jxbr) .and. va2max_jxb>0) lpencil_in(i_va2)=.true.
if (lpencil_in(i_jxbr)) then
lpencil_in(i_jxb)=.true.
lpencil_in(i_rho1)=.true.
endif
if (lpencil_in(i_jxb)) then
lpencil_in(i_jj)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_uxb2)) lpencil_in(i_uxb)=.true.
if (lpencil_in(i_uxb)) then
lpencil_in(i_uu)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_cosub)) then
lpencil_in(i_ub)=.true.
lpencil_in(i_u2)=.true.
lpencil_in(i_b2)=.true.
endif
if (lpencil_in(i_ub)) then
lpencil_in(i_uu)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_beta)) then
lpencil_in(i_b2)=.true.
lpencil_in(i_pp)=.true.
endif
if (lpencil_in(i_b2)) lpencil_in(i_bb)=.true.
if (lpencil_in(i_jj)) lpencil_in(i_bij)=.true.
if (lpencil_in(i_bb)) then
lpencil_in(i_aij)=.true.
endif
if (lpencil_in(i_djuidjbi)) then
lpencil_in(i_uij)=.true.
lpencil_in(i_bij)=.true.
endif
if (lpencil_in(i_jo)) then
lpencil_in(i_oo)=.true.
lpencil_in(i_jj)=.true.
endif
if (lpencil_in(i_ujxb)) then
lpencil_in(i_uu)=.true.
lpencil_in(i_jxb)=.true.
endif
if (lpencil_in(i_oxu)) then
lpencil_in(i_oo)=.true.
lpencil_in(i_uu)=.true.
endif
if (lpencil_in(i_oxuxb)) then
lpencil_in(i_oxu)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_jxbxb)) then
lpencil_in(i_jxb)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_jxbrxb)) then
lpencil_in(i_jxbr)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_glnrhoxb)) then
lpencil_in(i_glnrho)=.true.
lpencil_in(i_bb)=.true.
endif
if (lpencil_in(i_oxj)) then
lpencil_in(i_oo)=.true.
lpencil_in(i_jj)=.true.
endif
if (lpencil_in(i_jij)) lpencil_in(i_bij)=.true.
if (lpencil_in(i_sj)) then
lpencil_in(i_sij)=.true.
lpencil_in(i_jij)=.true.
endif
if (lpencil_in(i_ss12)) lpencil_in(i_sij)=.true.
if (lpencil_in(i_del2A)) then
lpencil_in(i_jj)=.true.
lpencil_in(i_graddivA)=.true.
endif
if (lpencil_in(i_ss12)) lpencil_in(i_sj)=.true.
!
endsubroutine pencil_interdep_magnetic
!***********************************************************************
subroutine calc_pencils_magnetic(f,p)
!
! Calculate Magnetic pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 19-nov-04/anders: coded
!
use Sub
use Deriv
use Diagnostics, only: sum_mn_name
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
! real, dimension (nx,3) :: bb_ext_pot
real, dimension (nx) :: rho1_jxb
real, dimension (nx) :: jcrossb2
real :: B2_ext,c,s
integer :: i,j,ix
!
intent(inout) :: f,p
! aa
if (lpencil(i_aa)) p%aa=f(l1:l2,m,n,iax:iaz)
! a2
if (lpencil(i_a2)) call dot2_mn(p%aa,p%a2)
! aij
if (lpencil(i_aij)) call gij(f,iaa,p%aij,1)
! diva
if (lpencil(i_diva)) call div_mn(p%aij,p%diva)
! bb
if (lpencil(i_bb)) then
call curl_mn(p%aij,p%bb)
!
! Save field before adding imposed field (for diagnostics).
!
p%bbb=p%bb
B2_ext=B_ext(1)**2+B_ext(2)**2+B_ext(3)**2
!
! Allow external field to precess about z-axis with frequency omega_Bz_ext.
!
if (B2_ext/=0.0) then
!
! luse_curvilinear_bext is default. The B_ext the user defines in
! magnetic_init_pars respects the coordinate system of preference
! which means that B_ext=(0.0,1.0,0.0) is an azimuthal field in cylindrical
! coordinates and a polar one in spherical.
!
if (omega_Bz_ext==0.0) then
B_ext_tmp=B_ext
elseif (omega_Bz_ext/=0.0) then
c=cos(omega_Bz_ext*t)
s=sin(omega_Bz_ext*t)
B_ext_tmp(1)=B_ext(1)*c-B_ext(2)*s
B_ext_tmp(2)=B_ext(1)*s+B_ext(2)*c
B_ext_tmp(3)=B_ext(3)
endif
!
! Add the external field.
!
if (B_ext_tmp(1)/=0.0) p%bb(:,1)=p%bb(:,1)+B_ext_tmp(1)
if (B_ext_tmp(2)/=0.0) p%bb(:,2)=p%bb(:,2)+B_ext_tmp(2)
if (B_ext_tmp(3)/=0.0) p%bb(:,3)=p%bb(:,3)+B_ext_tmp(3)
if (headtt) print*,'calc_pencils_magnetic: B_ext=',B_ext
if (headtt) print*,'calc_pencils_magnetic: B_ext_tmp=',B_ext_tmp
endif
!
! Add the external potential field.
!
if (lB_ext_pot) then
! call get_global(bb_ext_pot,m,n,'B_ext_pot')
! p%bb=p%bb+bb_ext_pot
endif
!
! Add external B-field.
!
if (iglobal_bx_ext/=0) p%bb(:,1)=p%bb(:,1)+f(l1:l2,m,n,iglobal_bx_ext)
if (iglobal_by_ext/=0) p%bb(:,2)=p%bb(:,2)+f(l1:l2,m,n,iglobal_by_ext)
if (iglobal_bz_ext/=0) p%bb(:,3)=p%bb(:,3)+f(l1:l2,m,n,iglobal_bz_ext)
endif
!
! b2 (default is that B_ext is not included), but this can be changed
! by setting luse_Bext_in_b2=.true.
!
if (luse_Bext_in_b2) then
if (lpencil(i_b2)) call dot2_mn(p%bb,p%b2)
else
if (lpencil(i_b2)) call dot2_mn(p%bbb,p%b2)
endif
! ab
if (lpencil(i_ab)) call dot_mn(p%aa,p%bbb,p%ab)
if (lpencil(i_ua)) call dot_mn(p%uu,p%aa,p%ua)
! uxb
if (lpencil(i_uxb)) then
call cross_mn(p%uu,p%bb,p%uxb)
call cross_mn(p%uu,p%bbb,uxbb)
! add external e-field.
if (iglobal_ex_ext/=0) p%uxb(:,1)=p%uxb(:,1)+f(l1:l2,m,n,iglobal_ex_ext)
if (iglobal_ey_ext/=0) p%uxb(:,2)=p%uxb(:,2)+f(l1:l2,m,n,iglobal_ey_ext)
if (iglobal_ez_ext/=0) p%uxb(:,3)=p%uxb(:,3)+f(l1:l2,m,n,iglobal_ez_ext)
endif
! exa
if (lpencil(i_exa)) call cross_mn(-p%uxb+eta*p%jj,p%aa,p%exa)
!
! bij, del2a, graddiva
! For non-cartesian coordinates jj is always required for del2a=graddiva-jj
!
if (lpencil(i_bij) .or. lpencil(i_del2a) .or. lpencil(i_graddiva) .or. &
lpencil(i_jj) ) then
call gij_etc(f,iaa,p%bij,p%del2a,p%graddiva)
if (.not. lpencil(i_bij)) p%bij=0.0 ! Avoid warnings from pencil
if (.not. lpencil(i_del2A)) p%del2A=0.0 ! consistency check...
if (.not. lpencil(i_graddiva)) p%graddiva=0.0
if (lpencil(i_jj)) call curl_mn(p%bij,p%jj)
endif
! jj
if (lpencil(i_jj)) then
p%jj=mu01*p%jj
!
! Add external j-field.
!
if (iglobal_jx_ext/=0) p%jj(:,1)=p%jj(:,1)+f(l1:l2,m,n,iglobal_jx_ext)
if (iglobal_jy_ext/=0) p%jj(:,2)=p%jj(:,2)+f(l1:l2,m,n,iglobal_jy_ext)
if (iglobal_jz_ext/=0) p%jj(:,3)=p%jj(:,3)+f(l1:l2,m,n,iglobal_jz_ext)
endif
! j2
if (lpencil(i_j2)) call dot2_mn(p%jj,p%j2)
! jb
if (lpencil(i_jb)) call dot_mn(p%jj,p%bbb,p%jb)
! va2
if (lpencil(i_va2)) then
p%va2=p%b2*mu01*p%rho1
if (lcheck_positive_va2 .and. minval(p%va2)<0.0) then
print*, 'calc_pencils_magnetic: Alfven speed is imaginary!'
print*, 'calc_pencils_magnetic: it, itsub, iproc=', it, itsub, iproc
print*, 'calc_pencils_magnetic: m, y(m), n, z(n)=', m, y(m), n, z(n)
p%va2=abs(p%va2)
endif
endif
! eta_va
if (lpencil(i_etava)) then
p%etava = mu0 * eta_va * dxmax * sqrt(p%va2)
if (eta_min > 0.) where (p%etava < eta_min) p%etava = 0.
endif
! eta_j
if (lpencil(i_etaj)) then
p%etaj = mu0 * eta_j * dxmax**2 * sqrt(mu0 * p%j2 * p%rho1)
if (eta_min > 0.) where (p%etaj < eta_min) p%etaj = 0.
endif
! eta_j2
if (lpencil(i_etaj2)) then
p%etaj2 = etaj20 * p%j2 * p%rho1
if (eta_min > 0.) where (p%etaj2 < eta_min) p%etaj2 = 0.
endif
! eta_jrho
if (lpencil(i_etajrho)) then
p%etajrho = mu0 * eta_jrho * dxmax * sqrt(p%j2) * p%rho1
if (eta_min > 0.) where (p%etajrho < eta_min) p%etajrho = 0.
endif
! jxb
if (lpencil(i_jxb)) call cross_mn(p%jj,p%bb,p%jxb)
! jxbr
if (lpencil(i_jxbr)) rho1_jxb=p%rho1
! cosjb
if (lpencil(i_cosjb)) then
do ix=1,nx
if ((abs(p%j2(ix)).le.tini).or.(abs(p%b2(ix)).le.tini))then
p%cosjb(ix)=0.
else
p%cosjb(ix)=p%jb(ix)/sqrt(p%j2(ix)*p%b2(ix))
endif
enddo
if (lpencil_check_at_work) then
! map penc0 value back to interval [-1,1]
p%cosjb = modulo(p%cosjb + 1.0, 2.0) - 1
endif
endif
! jparallel and jperp
if (lpencil(i_jparallel).or.lpencil(i_jperp)) then
p%jparallel=sqrt(p%j2)*p%cosjb
call dot2_mn(p%jxb,jcrossb2)
do ix=1,nx
if ((abs(p%j2(ix)).le.tini).or.(abs(p%b2(ix)).le.tini))then
p%jperp=0
else
p%jperp=sqrt(jcrossb2(ix))/sqrt(p%b2(ix))
endif
enddo
endif
! jxbr
if (lpencil(i_jxbr)) then
rho1_jxb=p%rho1
!
! Set rhomin_jxb>0 in order to limit the jxb term at very low densities.
! Set va2max_jxb>0 in order to limit the jxb term at very high Alfven speeds.
! Set va2power_jxb to an integer value in order to specify the power of the
! limiting term,
!
if (rhomin_jxb>0) rho1_jxb=min(rho1_jxb,1/rhomin_jxb)
if (va2max_jxb>0) then
rho1_jxb = rho1_jxb &
* (1+(p%va2/va2max_jxb)**va2power_jxb)**(-1.0/va2power_jxb)
endif
call multsv_mn(rho1_jxb,p%jxb,p%jxbr)
endif
! jxbr2
if (lpencil(i_jxbr2)) call dot2_mn(p%jxbr,p%jxbr2)
! ub
if (lpencil(i_ub)) call dot_mn(p%uu,p%bb,p%ub)
! cosub
if (lpencil(i_cosub)) then
do ix=1,nx
if ((abs(p%u2(ix)).le.tini).or.(abs(p%b2(ix)).le.tini)) then
p%cosub(ix)=0.
else
p%cosub(ix)=p%ub(ix)/sqrt(p%u2(ix)*p%b2(ix))
endif
enddo
if (lpencil_check) then
! map penc0 value back to interval [-1,1]
p%cosub = modulo(p%cosub + 1.0, 2.0) - 1
endif
endif
! uxb2
if (lpencil(i_uxb2)) call dot2_mn(p%uxb,p%uxb2)
! uxj
if (lpencil(i_uxj)) call cross_mn(p%uu,p%jj,p%uxj)
! beta
if (lpencil(i_beta)) p%beta=0.5*p%b2/p%pp
! djuidjbi
if (lpencil(i_djuidjbi)) call multmm_sc(p%uij,p%bij,p%djuidjbi)
! jo
if (lpencil(i_jo)) call dot(p%jj,p%oo,p%jo)
! ujxb
if (lpencil(i_ujxb)) call dot_mn(p%uu,p%jxb,p%ujxb)
! oxu
if (lpencil(i_oxu)) call cross_mn(p%oo,p%uu,p%oxu)
! oxuxb
if (lpencil(i_oxuxb)) call cross_mn(p%oxu,p%bb,p%oxuxb)
! jxbxb
if (lpencil(i_jxbxb)) call cross_mn(p%jxb,p%bb,p%jxbxb)
! jxbrxb
if (lpencil(i_jxbrxb)) call cross_mn(p%jxbr,p%bb,p%jxbrxb)
! glnrhoxb
if (lpencil(i_glnrhoxb)) call cross_mn(p%glnrho,p%bb,p%glnrhoxb)
! oxj
if (lpencil(i_oxj)) call cross_mn(p%oo,p%jj,p%oxJ)
! jij
if (lpencil(i_jij)) then
do j=1,3
do i=1,3
p%jij(:,i,j)=.5*(p%bij(:,i,j)+p%bij(:,j,i))
enddo
enddo
endif
! sj
if (lpencil(i_sj)) call multmm_sc(p%sij,p%jij,p%sj)
! ss12
if (lpencil(i_ss12)) p%ss12=sqrt(abs(p%sj))
!
! Store bb in auxiliary variable if requested.
! Just neccessary immediately before writing snapshots, but how would we
! know we are?
!
if (lbb_as_aux) f(l1:l2,m,n,ibx:ibz)=p%bb
if (ljj_as_aux) f(l1:l2,m,n,ijx:ijz)=p%jj
!
endsubroutine calc_pencils_magnetic
!***********************************************************************
subroutine daa_dt(df,p)
!
! Magnetic field evolution.
!
! Calculate dA/dt=uxB+3/2 Omega_0 A_y x_dir -eta mu_0 J.
! Add jxb/rho to momentum equation.
! Add eta mu_0 j2/rho to entropy equation.
!
! 22-nov-01/nils: coded
! 1-may-02/wolf: adapted for pencil_modular
! 17-jun-03/ulf: added bx^2, by^2 and bz^2 as separate diagnostics
! 8-aug-03/axel: introduced B_ext21=1./B_ext**2, and set =1 to avoid div. by 0
! 26-may-04/axel: ambipolar diffusion added
!
use Diagnostics
use EquationOfState, only: eoscalc,gamma_m1,cs0
use Io, only: output_pencil
use Mpicomm, only: stop_it
use Sub
!
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
real, dimension (nx,3) :: uxDxuxb,fres
real, dimension (nx,3) :: exj,dexb,phib,jxbb
real, dimension (nx) :: exabot,exatop
real, dimension (nx) :: jxb_dotB0,uxb_dotB0
real, dimension (nx) :: oxuxb_dotB0,jxbxb_dotB0,uxDxuxb_dotB0
real, dimension (nx) :: uj,aj,phi
real, dimension (nx) :: uxj_dotB0,b3b21,b1b32,b2b13
real, dimension (nx) :: rho1_jxb
real, dimension (nx) :: B1dot_glnrhoxb
real, dimension (nx) :: eta_smag,etatotal
real, dimension (nx) :: fres2
real, parameter :: OmegaSS=1.0
integer :: i
integer, parameter :: nxy=nxgrid*nygrid
!
intent(inout) :: p
intent(inout) :: df
!
! Identify module and boundary conditions.
!
call timing('daa_dt','entered',mnloop=.true.)
if (headtt.or.ldebug) print*,'daa_dt: SOLVE'
if (headtt) then
call identify_bcs('Ax',iax)
call identify_bcs('Ay',iay)
call identify_bcs('Az',iaz)
endif
!
! Add jxb/rho to momentum equation.
!
if (lhydro.and.llorentzforce) &
df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+p%jxbr
!
! Restivivity term
!
! Because of gauge invariance, we can add the gradient of an arbitrary scalar
! field Phi to the induction equation without changing the magnetic field,
! dA/dt = u x B - eta j + grad(Phi).
!
! If lweyl_gauge=T, we choose Phi = const. and solve
! Else, if lweyl_gauge=F, we make the gauge choice Phi = eta div(A)
! and thus solve
! dA/dt = u x B + eta laplace(A) + div(A) grad(eta).
!
! Note: lweyl_gauge=T is so far only implemented for some resistivity types.
!
if (headtt) print*, 'daa_dt: iresistivity=', iresistivity
!
fres=0.0
etatotal=0.0
!
if (lresi_eta_const) then
fres=fres+eta*p%del2a
if (lfirst.and.ldt) diffus_eta=diffus_eta+eta
etatotal=etatotal+eta
endif
!
if (lresi_eta_shock) then
do i=1,3
fres(:,i)=fres(:,i)+ &
eta_shock*(p%shock*p%del2a(:,i)+p%diva*p%gshock(:,i))
enddo
if (lfirst.and.ldt) diffus_eta=diffus_eta+eta_shock*p%shock
etatotal=etatotal+eta_shock*p%shock
endif
!
! Multiply resistivity by Nyquist scale, for resistive time-step.
!
if (lfirst.and.ldt) then
diffus_eta =diffus_eta*dxyz_2
if (headtt.or.ldebug) then
print*, 'daa_dt: max(diffus_eta) =', maxval(diffus_eta)
endif
endif
!
! Add Ohmic heat to entropy or temperature equation.
!
if (lohmic_heat.and.lentropy) &
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + &
etatotal*mu0*p%j2*p%rho1*p%TT1
!
! Induction equation.
!
if (linduction) &
df(l1:l2,m,n,iax:iaz)=df(l1:l2,m,n,iax:iaz)+p%uxb+fres
!
! Add ``va^2/dx^2'' contribution to timestep.
! Consider advective timestep only when lhydro=T.
!
if (lfirst.and.ldt) then
if (lhydro) then
rho1_jxb=p%rho1
if (rhomin_jxb>0) rho1_jxb=min(rho1_jxb,1/rhomin_jxb)
if (va2max_jxb>0) then
rho1_jxb = rho1_jxb &
* (1+(p%va2/va2max_jxb)**va2power_jxb)**(-1.0/va2power_jxb)
endif
advec_va2=((p%bb(:,1)*dx_1(l1:l2))**2+ &
(p%bb(:,2)*dy_1( m ))**2+ &
(p%bb(:,3)*dz_1( n ))**2)*mu01*rho1_jxb
endif
endif
!
! Calculate diagnostic quantities.
!
if (ldiagnos) then
if (idiag_beta1m/=0) call sum_mn_name(p%beta,idiag_beta1m)
if (idiag_beta1max/=0) call max_mn_name(p%beta,idiag_beta1max)
!
! Contributions to vertical Poynting vector. Consider them here
! separately, because the contribution b^2*uz may be a good indicator
! of magnetic buoyancy.
!
if (idiag_b2ruzm/=0) call sum_mn_name(p%b2*p%rho*p%uu(:,3),idiag_b2ruzm)
if (idiag_b2uzm/=0) call sum_mn_name(p%b2*p%uu(:,3),idiag_b2uzm)
if (idiag_ubbzm/=0) call sum_mn_name(p%ub*p%bb(:,3),idiag_ubbzm)
!
! Mean squared and maximum squared magnetic field.
!
if (idiag_b1m/=0) call sum_mn_name(sqrt(p%b2),idiag_b1m)
if (idiag_b2m/=0) call sum_mn_name(p%b2,idiag_b2m)
if (idiag_bm2/=0) call max_mn_name(p%b2,idiag_bm2)
if (idiag_brms/=0) call sum_mn_name(p%b2,idiag_brms,lsqrt=.true.)
if (idiag_emag/=0) call integrate_mn_name(mu012*p%b2,idiag_emag)
if (idiag_brmsh/=0) then
if (lequatory) call sum_mn_name_halfy(p%b2,idiag_brmsh)
if (lequatorz) call sum_mn_name_halfz(p%b2,idiag_brmsh)
fname(idiag_brmsn)=fname_half(idiag_brmsh,1)
fname(idiag_brmss)=fname_half(idiag_brmsh,2)
itype_name(idiag_brmsn)=ilabel_sum_sqrt
itype_name(idiag_brmss)=ilabel_sum_sqrt
endif
if (idiag_bmax/=0) call max_mn_name(p%b2,idiag_bmax,lsqrt=.true.)
if (idiag_bxmin/=0) call max_mn_name(-p%bb(:,1),idiag_bxmin,lneg=.true.)
if (idiag_bymin/=0) call max_mn_name(-p%bb(:,2),idiag_bymin,lneg=.true.)
if (idiag_bzmin/=0) call max_mn_name(-p%bb(:,3),idiag_bzmin,lneg=.true.)
if (idiag_bxmax/=0) call max_mn_name(p%bb(:,1),idiag_bxmax)
if (idiag_bymax/=0) call max_mn_name(p%bb(:,2),idiag_bymax)
if (idiag_bzmax/=0) call max_mn_name(p%bb(:,3),idiag_bzmax)
if (idiag_aybym2/=0) &
call sum_mn_name(2*p%aa(:,2)*p%bb(:,2),idiag_aybym2)
if (idiag_abm/=0) call sum_mn_name(p%ab,idiag_abm)
if (idiag_abumx/=0) call sum_mn_name(p%uu(:,1)*p%ab,idiag_abumx)
if (idiag_abumy/=0) call sum_mn_name(p%uu(:,2)*p%ab,idiag_abumy)
if (idiag_abumz/=0) call sum_mn_name(p%uu(:,3)*p%ab,idiag_abumz)
if (idiag_abrms/=0) call sum_mn_name(p%ab**2,idiag_abrms,lsqrt=.true.)
!
! Hemispheric magnetic helicity of total field.
! North means 1 and south means 2.
!
if (idiag_abmh/=0) then
if (lequatory) call sum_mn_name_halfy(p%ab,idiag_abmh)
if (lequatorz) call sum_mn_name_halfz(p%ab,idiag_abmh)
fname(idiag_abmn)=fname_half(idiag_abmh,1)
fname(idiag_abms)=fname_half(idiag_abmh,2)
itype_name(idiag_abmn)=ilabel_sum
itype_name(idiag_abms)=ilabel_sum
endif
!
! Hemispheric current helicity of total field.
! North means 1 and south means 2.
!
if (idiag_jbmh/=0) then
if (lequatory) call sum_mn_name_halfy(p%jb,idiag_jbmh)
if (lequatorz) call sum_mn_name_halfz(p%jb,idiag_jbmh)
fname(idiag_jbmn)=fname_half(idiag_jbmh,1)
fname(idiag_jbms)=fname_half(idiag_jbmh,2)
itype_name(idiag_jbmn)=ilabel_sum
itype_name(idiag_jbms)=ilabel_sum
endif
!
! Cross helicity (linkage between vortex tubes and flux tubes).
!
if (idiag_ubm/=0) call sum_mn_name(p%ub,idiag_ubm)
if (idiag_uxbxm/=0) call sum_mn_name(p%uu(:,1)*p%bb(:,1),idiag_uxbxm)
if (idiag_uybym/=0) call sum_mn_name(p%uu(:,2)*p%bb(:,2),idiag_uybym)
if (idiag_uzbzm/=0) call sum_mn_name(p%uu(:,3)*p%bb(:,3),idiag_uzbzm)
if (idiag_cosubm/=0) call sum_mn_name(p%cosub,idiag_cosubm)
!
! Field-velocity cross helicity (linkage between velocity and magnetic tubes).
!
if (idiag_uam/=0) call sum_mn_name(p%ua,idiag_uam)
!
! Current-vortex cross helicity (linkage between vortex and current tubes).
!
if (idiag_ujm/=0) then
call dot (p%uu,p%jj,uj)
call sum_mn_name(uj,idiag_ujm)
endif
!
! rms values of mean field.
!
if (idiag_bxm/=0) call sum_mn_name(p%bbb(:,1),idiag_bxm)
if (idiag_bym/=0) call sum_mn_name(p%bbb(:,2),idiag_bym)
if (idiag_bzm/=0) call sum_mn_name(p%bbb(:,3),idiag_bzm)
if (idiag_bx2m/=0) call sum_mn_name(p%bbb(:,1)**2,idiag_bx2m)
if (idiag_by2m/=0) call sum_mn_name(p%bbb(:,2)**2,idiag_by2m)
if (idiag_bz2m/=0) call sum_mn_name(p%bbb(:,3)**2,idiag_bz2m)
if (idiag_bxbym/=0) call sum_mn_name(p%bbb(:,1)*p%bbb(:,2),idiag_bxbym)
if (idiag_bxbzm/=0) call sum_mn_name(p%bbb(:,1)*p%bbb(:,3),idiag_bxbzm)
if (idiag_bybzm/=0) call sum_mn_name(p%bbb(:,2)*p%bbb(:,3),idiag_bybzm)
if (idiag_djuidjbim/=0) call sum_mn_name(p%djuidjbi,idiag_djuidjbim)
!
! Calculate B*sin(phi) = -<Bx*sinkz> + <By*coskz>.
!
if (idiag_bsinphz/=0) call sum_mn_name(-p%bbb(:,1)*sinkz(n)+p%bbb(:,2)*coskz(n),idiag_bsinphz)
!
! Calculate B*cos(phi) = <Bx*coskz> + <By*sinkz>.
!
if (idiag_bcosphz/=0) call sum_mn_name(+p%bbb(:,1)*coskz(n)+p%bbb(:,2)*sinkz(n),idiag_bcosphz)
!
! Magnetic field components at one point (=pt).
!
if (lroot.and.m==mpoint.and.n==npoint) then
if (idiag_bxpt/=0) call save_name(p%bb(lpoint-nghost,1),idiag_bxpt)
if (idiag_bypt/=0) call save_name(p%bb(lpoint-nghost,2),idiag_bypt)
if (idiag_bzpt/=0) call save_name(p%bb(lpoint-nghost,3),idiag_bzpt)
if (idiag_Expt/=0) call save_name(uxbb(lpoint-nghost,1),idiag_Expt)
if (idiag_Eypt/=0) call save_name(uxbb(lpoint-nghost,2),idiag_Eypt)
if (idiag_Ezpt/=0) call save_name(uxbb(lpoint-nghost,3),idiag_Ezpt)
endif
!
! v_A = |B|/sqrt(rho); in units where mu_0=1
!
if (idiag_vA2m/=0) call sum_mn_name(p%va2,idiag_vA2m)
if (idiag_vArms/=0) call sum_mn_name(p%va2,idiag_vArms,lsqrt=.true.)
if (idiag_vAmax/=0) call max_mn_name(p%va2,idiag_vAmax,lsqrt=.true.)
if (idiag_dtb/=0) &
call max_mn_name(sqrt(advec_va2)/cdt,idiag_dtb,l_dt=.true.)
!
! Lorentz force.
!
if (idiag_jxbrxm/=0) call sum_mn_name(p%jxbr(:,1),idiag_jxbrxm)
if (idiag_jxbrym/=0) call sum_mn_name(p%jxbr(:,2),idiag_jxbrym)
if (idiag_jxbrzm/=0) call sum_mn_name(p%jxbr(:,3),idiag_jxbrzm)
if (idiag_jxbr2m/=0) call sum_mn_name(p%jxbr2,idiag_jxbr2m)
!
! <J.A> for calculating k_effective, for example.
!
if (idiag_ajm/=0) then
call dot (p%aa,p%jj,aj)
call sum_mn_name(aj,idiag_ajm)
endif
!
! Output kx_aa for calculating k_effective.
!
if (idiag_kx_aa/=0) call save_name(kx_aa(1),idiag_kx_aa)
!
! Helicity integrals.
!
if (idiag_ab_int/=0) call integrate_mn_name(p%ab,idiag_ab_int)
if (idiag_jb_int/=0) call integrate_mn_name(p%jb,idiag_jb_int)
!
! <J.B>
!
if (idiag_jbm/=0) call sum_mn_name(p%jb,idiag_jbm)
if (idiag_j2m/=0) call sum_mn_name(p%j2,idiag_j2m)
if (idiag_jm2/=0) call max_mn_name(p%j2,idiag_jm2)
if (idiag_jrms/=0) call sum_mn_name(p%j2,idiag_jrms,lsqrt=.true.)
if (idiag_jmax/=0) call max_mn_name(p%j2,idiag_jmax,lsqrt=.true.)
if (idiag_epsM_LES/=0) call sum_mn_name(eta_smag*p%j2,idiag_epsM_LES)
if (idiag_dteta/=0) call max_mn_name(diffus_eta/cdtv,idiag_dteta,l_dt=.true.)
if (idiag_cosjbm/=0) call sum_mn_name(p%cosjb,idiag_cosjbm)
if (idiag_jparallelm/=0) call sum_mn_name(p%jparallel,idiag_jparallelm)
if (idiag_jperpm/=0) call sum_mn_name(p%jperp,idiag_jperpm)
!
! Resistivity.
!
if (idiag_etasmagm/=0) call sum_mn_name(eta_smag,idiag_etasmagm)
if (idiag_etasmagmin/=0) call max_mn_name(-eta_smag,idiag_etasmagmin,lneg=.true.)
if (idiag_etasmagmax/=0) call max_mn_name(eta_smag,idiag_etasmagmax)
if (idiag_etavamax/=0) call max_mn_name(p%etava,idiag_etavamax)
if (idiag_etajmax/=0) call max_mn_name(p%etaj,idiag_etajmax)
if (idiag_etaj2max/=0) call max_mn_name(p%etaj2,idiag_etaj2max)
if (idiag_etajrhomax/=0) call max_mn_name(p%etajrho,idiag_etajrhomax)
!
! Not correct for hyperresistivity:
!
if (idiag_epsM/=0) call sum_mn_name(eta*p%j2,idiag_epsM)
!
! Heating by ion-neutrals friction.
!
if (idiag_epsAD/=0) call sum_mn_name(nu_ni1*p%rho*p%jxbr2,idiag_epsAD)
!
! <A>'s, <A^2> and A^2|max
!
if (idiag_axm/=0) call sum_mn_name(p%aa(:,1),idiag_axm)
if (idiag_aym/=0) call sum_mn_name(p%aa(:,2),idiag_aym)
if (idiag_azm/=0) call sum_mn_name(p%aa(:,3),idiag_azm)
if (idiag_arms/=0) call sum_mn_name(p%a2,idiag_arms,lsqrt=.true.)
if (idiag_amax/=0) call max_mn_name(p%a2,idiag_amax,lsqrt=.true.)
!
if (idiag_uxbm /=0 .or. &
idiag_uxbmx/=0 .or. idiag_uxbmy/=0 .or. idiag_uxbmz/=0) then
call dot(B_ext_inv,p%uxb,uxb_dotB0)
if (idiag_uxbm/=0) call sum_mn_name(uxb_dotB0,idiag_uxbm)
if (idiag_uxbmx/=0) call sum_mn_name(uxbb(:,1),idiag_uxbmx)
if (idiag_uxbmy/=0) call sum_mn_name(uxbb(:,2),idiag_uxbmy)
if (idiag_uxbmz/=0) call sum_mn_name(uxbb(:,3),idiag_uxbmz)
endif
!
! Calculate part I of magnetic helicity flux (ExA contribution).
!
if (idiag_examx/=0 .or. idiag_examy/=0 .or. idiag_examz/=0 .or. &
idiag_exatop/=0 .or. idiag_exabot/=0) then
if (idiag_examx/=0) call sum_mn_name(p%exa(:,1),idiag_examx)
if (idiag_examy/=0) call sum_mn_name(p%exa(:,2),idiag_examy)
if (idiag_examz/=0) call sum_mn_name(p%exa(:,3),idiag_examz)
!
if (idiag_exabot/=0) then
if (z(n)==xyz0(3)) then
exabot=p%exa(:,3)
else
exabot=0.
endif
call integrate_mn_name(exabot,idiag_exabot)
endif
!
if (idiag_exatop/=0) then
if (z(n)==xyz1(3)) then
exatop=p%exa(:,3)
else
exatop=0.
endif
call integrate_mn_name(exatop,idiag_exatop)
endif
!
endif
!
! Calculate part II of magnetic helicity flux (phi*B contribution).
!
if (idiag_phibmx/=0 .or. idiag_phibmy/=0 .or. idiag_phibmz/=0) then
phi=eta*p%diva
call multvs(p%bb,phi,phib)
if (idiag_phibmx/=0) call sum_mn_name(phib(:,1),idiag_phibmx)
if (idiag_phibmy/=0) call sum_mn_name(phib(:,2),idiag_phibmy)
if (idiag_phibmz/=0) call sum_mn_name(phib(:,3),idiag_phibmz)
endif
!
! Calculate part I of current helicity flux (for imposed field).
!
if (idiag_exjmx/=0 .or. idiag_exjmy/=0 .or. idiag_exjmz/=0) then
call cross_mn(-p%uxb+eta*p%jj,p%jj,exj)
if (idiag_exjmx/=0) call sum_mn_name(exj(:,1),idiag_exjmx)
if (idiag_exjmy/=0) call sum_mn_name(exj(:,2),idiag_exjmy)
if (idiag_exjmz/=0) call sum_mn_name(exj(:,3),idiag_exjmz)
endif
!
! Calculate part II of current helicity flux (for imposed field).
! < curlE x B >|_i = < B_{j,i} E_j >
! Use the full B (with B_ext)
!
if (idiag_dexbmx/=0 .or. idiag_dexbmy/=0 .or. idiag_dexbmz/=0) then
call multmv_transp(p%bij,-p%uxb+eta*p%jj,dexb)
if (idiag_dexbmx/=0) call sum_mn_name(dexb(:,1),idiag_dexbmx)
if (idiag_dexbmy/=0) call sum_mn_name(dexb(:,2),idiag_dexbmy)
if (idiag_dexbmz/=0) call sum_mn_name(dexb(:,3),idiag_dexbmz)
endif
!
! Calculate <uxj>.B0/B0^2.
!
if (idiag_uxjm/=0) then
call dot(B_ext_inv,p%uxj,uxj_dotB0)
call sum_mn_name(uxj_dotB0,idiag_uxjm)
endif
!
! Calculate <u x B>_rms, <resistive terms>_rms, <ratio ~ Rm>_rms.
!
if (idiag_uxBrms/=0) call sum_mn_name(p%uxb2,idiag_uxBrms,lsqrt=.true.)
if (idiag_Bresrms/=0 .or. idiag_Rmrms/=0) then
call dot2_mn(fres,fres2)
if (idiag_Bresrms/=0) &
call sum_mn_name(fres2,idiag_Bresrms,lsqrt=.true.)
if (idiag_Rmrms/=0) &
call sum_mn_name(p%uxb2/fres2,idiag_Rmrms,lsqrt=.true.)
endif
!
! Calculate <b^2*divu>, which is part of <u.(jxb)>.
! Note that <u.(jxb)>=1/2*<b^2*divu>+<u.bgradb>.
!
if (idiag_b2divum/=0) call sum_mn_name(p%b2*p%divu,idiag_b2divum)
!
! Calculate <u.(jxb)>.
!
if (idiag_ujxbm/=0) call sum_mn_name(p%ujxb,idiag_ujxbm)
!
! Calculate <jxb>.B_0/B_0^2.
!
call dot(B_ext_inv,p%jxb,jxb_dotB0)
if (idiag_jxbm/=0) call sum_mn_name(jxb_dotB0,idiag_jxbm)
if (idiag_jxbmx/=0.or.idiag_jxbmy/=0.or.idiag_jxbmz/=0) then
call cross_mn(p%jj,p%bbb,jxbb)
if (idiag_jxbmx/=0) call sum_mn_name(jxbb(:,1),idiag_jxbmx)
if (idiag_jxbmy/=0) call sum_mn_name(jxbb(:,2),idiag_jxbmy)
if (idiag_jxbmz/=0) call sum_mn_name(jxbb(:,3),idiag_jxbmz)
endif
!
! Magnetic triple correlation term (for imposed field).
!
if (idiag_jxbxbm/=0) then
call dot(B_ext_inv,p%jxbxb,jxbxb_dotB0)
call sum_mn_name(jxbxb_dotB0,idiag_jxbxbm)
endif
!
! Triple correlation from Reynolds tensor (for imposed field).
!
if (idiag_oxuxbm/=0) then
call dot(B_ext_inv,p%oxuxb,oxuxb_dotB0)
call sum_mn_name(oxuxb_dotB0,idiag_oxuxbm)
endif
!
! Triple correlation from pressure gradient (for imposed field).
! (assume cs2=1, and that no entropy evolution is included)
! This is ok for all applications currently under consideration.
!
if (idiag_gpxbm/=0) then
call dot_mn_sv(B1_ext,p%glnrhoxb,B1dot_glnrhoxb)
call sum_mn_name(B1dot_glnrhoxb,idiag_gpxbm)
endif
!
! < u x curl(uxB) > = < E_i u_{j,j} - E_j u_{j,i} >
! ( < E_1 u2,2 + E1 u3,3 - E2 u2,1 - E3 u3,1 >
! < E_2 u1,1 + E2 u3,3 - E1 u2,1 - E3 u3,2 >
! < E_3 u1,1 + E3 u2,2 - E1 u3,1 - E2 u2,3 > )
!
if (idiag_uxDxuxbm/=0) then
uxDxuxb(:,1)=p%uxb(:,1)*(p%uij(:,2,2)+p%uij(:,3,3))-p%uxb(:,2)*p%uij(:,2,1)-p%uxb(:,3)*p%uij(:,3,1)
uxDxuxb(:,2)=p%uxb(:,2)*(p%uij(:,1,1)+p%uij(:,3,3))-p%uxb(:,1)*p%uij(:,1,2)-p%uxb(:,3)*p%uij(:,3,2)
uxDxuxb(:,3)=p%uxb(:,3)*(p%uij(:,1,1)+p%uij(:,2,2))-p%uxb(:,1)*p%uij(:,1,3)-p%uxb(:,2)*p%uij(:,2,3)
call dot(B_ext_inv,uxDxuxb,uxDxuxb_dotB0)
call sum_mn_name(uxDxuxb_dotB0,idiag_uxDxuxbm)
endif
!
! alpM11=<b3*b2,1>
!
if (idiag_b3b21m/=0) then
b3b21=p%bb(:,3)*p%bij(:,2,1)
call sum_mn_name(b3b21,idiag_b3b21m)
endif
!
! alpM22=<b1*b3,2>
!
if (idiag_b1b32m/=0) then
b1b32=p%bb(:,1)*p%bij(:,3,2)
call sum_mn_name(b1b32,idiag_b1b32m)
endif
!
! alpM33=<b2*b1,3>
!
if (idiag_b2b13m/=0) then
b2b13=p%bb(:,2)*p%bij(:,1,3)
call sum_mn_name(b2b13,idiag_b2b13m)
endif
!
endif ! endif (ldiagnos)
!
! 1d-averages. Happens at every it1d timesteps, NOT at every it1.
!
if (l1davgfirst .or. (ldiagnos .and. ldiagnos_need_zaverages)) then
if (idiag_bxmx/=0) call yzsum_mn_name_x(p%bb(:,1),idiag_bxmx)
if (idiag_bymx/=0) call yzsum_mn_name_x(p%bb(:,2),idiag_bymx)
if (idiag_bzmx/=0) call yzsum_mn_name_x(p%bb(:,3),idiag_bzmx)
if (idiag_bx2mx/=0) call yzsum_mn_name_x(p%bb(:,1)**2,idiag_bx2mx)
if (idiag_by2mx/=0) call yzsum_mn_name_x(p%bb(:,2)**2,idiag_by2mx)
if (idiag_bz2mx/=0) call yzsum_mn_name_x(p%bb(:,3)**2,idiag_bz2mx)
if (idiag_bxbymx/=0) &
call yzsum_mn_name_x(p%bbb(:,1)*p%bbb(:,2),idiag_bxbymx)
if (idiag_bxmy/=0) call xzsum_mn_name_y(p%bb(:,1),idiag_bxmy)
if (idiag_bymy/=0) call xzsum_mn_name_y(p%bb(:,2),idiag_bymy)
if (idiag_bzmy/=0) call xzsum_mn_name_y(p%bb(:,3),idiag_bzmy)
if (idiag_bx2my/=0) call xzsum_mn_name_y(p%bb(:,1)**2,idiag_bx2my)
if (idiag_by2my/=0) call xzsum_mn_name_y(p%bb(:,2)**2,idiag_by2my)
if (idiag_bz2my/=0) call xzsum_mn_name_y(p%bb(:,3)**2,idiag_bz2my)
if (idiag_axmz/=0) call xysum_mn_name_z(p%aa(:,1),idiag_axmz)
if (idiag_aymz/=0) call xysum_mn_name_z(p%aa(:,2),idiag_aymz)
if (idiag_azmz/=0) call xysum_mn_name_z(p%aa(:,3),idiag_azmz)
if (idiag_abuxmz/=0) call xysum_mn_name_z(p%ab*p%uu(:,1),idiag_abuxmz)
if (idiag_abuymz/=0) call xysum_mn_name_z(p%ab*p%uu(:,2),idiag_abuymz)
if (idiag_abuzmz/=0) call xysum_mn_name_z(p%ab*p%uu(:,3),idiag_abuzmz)
if (idiag_uabxmz/=0) call xysum_mn_name_z(p%ua*p%bb(:,1),idiag_uabxmz)
if (idiag_uabymz/=0) call xysum_mn_name_z(p%ua*p%bb(:,2),idiag_uabymz)
if (idiag_uabzmz/=0) call xysum_mn_name_z(p%ua*p%bb(:,3),idiag_uabzmz)
if (idiag_bxmz/=0) call xysum_mn_name_z(p%bb(:,1),idiag_bxmz)
if (idiag_bymz/=0) call xysum_mn_name_z(p%bb(:,2),idiag_bymz)
if (idiag_bzmz/=0) call xysum_mn_name_z(p%bb(:,3),idiag_bzmz)
if (idiag_jxmz/=0) call xysum_mn_name_z(p%jj(:,1),idiag_jxmz)
if (idiag_jymz/=0) call xysum_mn_name_z(p%jj(:,2),idiag_jymz)
if (idiag_jzmz/=0) call xysum_mn_name_z(p%jj(:,3),idiag_jzmz)
if (idiag_Exmz/=0) call xysum_mn_name_z(p%uxb(:,1),idiag_Exmz)
if (idiag_Eymz/=0) call xysum_mn_name_z(p%uxb(:,2),idiag_Eymz)
if (idiag_Ezmz/=0) call xysum_mn_name_z(p%uxb(:,3),idiag_Ezmz)
if (idiag_bx2mz/=0) call xysum_mn_name_z(p%bb(:,1)**2,idiag_bx2mz)
if (idiag_by2mz/=0) call xysum_mn_name_z(p%bb(:,2)**2,idiag_by2mz)
if (idiag_bz2mz/=0) call xysum_mn_name_z(p%bb(:,3)**2,idiag_bz2mz)
if (idiag_bx2rmz/=0) call xysum_mn_name_z(p%bb(:,1)**2*p%rho1,idiag_bx2rmz)
if (idiag_by2rmz/=0) call xysum_mn_name_z(p%bb(:,2)**2*p%rho1,idiag_by2rmz)
if (idiag_bz2rmz/=0) call xysum_mn_name_z(p%bb(:,3)**2*p%rho1,idiag_bz2rmz)
if (idiag_jbmz/=0) call xysum_mn_name_z(p%jb,idiag_jbmz)
if (idiag_abmz/=0) call xysum_mn_name_z(p%ab,idiag_abmz)
if (idiag_uamz/=0) call xysum_mn_name_z(p%ua,idiag_uamz)
if (idiag_etatotalmx/=0) call yzsum_mn_name_x(etatotal,idiag_etatotalmx)
if (idiag_etatotalmz/=0) call xysum_mn_name_z(etatotal,idiag_etatotalmz)
!
! Calculate magnetic helicity flux (ExA contribution).
!
if (idiag_examz1/=0) call xysum_mn_name_z(p%exa(:,1),idiag_examz1)
if (idiag_examz2/=0) call xysum_mn_name_z(p%exa(:,2),idiag_examz2)
if (idiag_examz3/=0) call xysum_mn_name_z(p%exa(:,3),idiag_examz3)
!
! Maxwell stress components.
!
if (idiag_bxbymy/=0) &
call xzsum_mn_name_y(p%bbb(:,1)*p%bbb(:,2),idiag_bxbymy)
if (idiag_bxbzmy/=0) &
call xzsum_mn_name_y(p%bbb(:,1)*p%bbb(:,3),idiag_bxbzmy)
if (idiag_bybzmy/=0) &
call xzsum_mn_name_y(p%bbb(:,2)*p%bbb(:,3),idiag_bybzmy)
if (idiag_bxbymz/=0) &
call xysum_mn_name_z(p%bbb(:,1)*p%bbb(:,2),idiag_bxbymz)
if (idiag_bxbzmz/=0) &
call xysum_mn_name_z(p%bbb(:,1)*p%bbb(:,3),idiag_bxbzmz)
if (idiag_bybzmz/=0) &
call xysum_mn_name_z(p%bbb(:,2)*p%bbb(:,3),idiag_bybzmz)
if (idiag_jxbrxmx/=0) call yzsum_mn_name_x(p%jxbr(:,1),idiag_jxbrxmx)
if (idiag_jxbrymx/=0) call yzsum_mn_name_x(p%jxbr(:,2),idiag_jxbrymx)
if (idiag_jxbrzmx/=0) call yzsum_mn_name_x(p%jxbr(:,3),idiag_jxbrzmx)
if (idiag_jxbrxmy/=0) call xzsum_mn_name_y(p%jxbr(:,1),idiag_jxbrxmy)
if (idiag_jxbrymy/=0) call xzsum_mn_name_y(p%jxbr(:,2),idiag_jxbrymy)
if (idiag_jxbrzmy/=0) call xzsum_mn_name_y(p%jxbr(:,3),idiag_jxbrzmy)
if (idiag_jxbrxmz/=0) call xysum_mn_name_z(p%jxbr(:,1),idiag_jxbrxmz)
if (idiag_jxbrymz/=0) call xysum_mn_name_z(p%jxbr(:,2),idiag_jxbrymz)
if (idiag_jxbrzmz/=0) call xysum_mn_name_z(p%jxbr(:,3),idiag_jxbrzmz)
if (idiag_b2mz/=0) call xysum_mn_name_z(p%b2,idiag_b2mz)
if (idiag_mflux_x/=0) &
call yzintegrate_mn_name_x(p%bb(:,1),idiag_mflux_x)
if (idiag_mflux_y/=0) &
call xzintegrate_mn_name_y(p%bb(:,2),idiag_mflux_y)
if (idiag_mflux_z/=0) &
call xyintegrate_mn_name_z(p%bb(:,3),idiag_mflux_z)
endif
!
! 2-D averages.
! Note that this does not necessarily happen with ldiagnos=.true.
!
if (l2davgfirst) then
if (idiag_bxmxy/=0) call zsum_mn_name_xy(p%bb(:,1),idiag_bxmxy)
if (idiag_bymxy/=0) call zsum_mn_name_xy(p%bb(:,2),idiag_bymxy)
if (idiag_bzmxy/=0) call zsum_mn_name_xy(p%bb(:,3),idiag_bzmxy)
if (idiag_jxmxy/=0) call zsum_mn_name_xy(p%jj(:,1),idiag_jxmxy)
if (idiag_jymxy/=0) call zsum_mn_name_xy(p%jj(:,2),idiag_jymxy)
if (idiag_jzmxy/=0) call zsum_mn_name_xy(p%jj(:,3),idiag_jzmxy)
if (idiag_axmxz/=0) call ysum_mn_name_xz(p%aa(:,1),idiag_axmxz)
if (idiag_aymxz/=0) call ysum_mn_name_xz(p%aa(:,2),idiag_aymxz)
if (idiag_azmxz/=0) call ysum_mn_name_xz(p%aa(:,3),idiag_azmxz)
if (idiag_bxmxz/=0) call ysum_mn_name_xz(p%bb(:,1),idiag_bxmxz)
if (idiag_bymxz/=0) call ysum_mn_name_xz(p%bb(:,2),idiag_bymxz)
if (idiag_bzmxz/=0) call ysum_mn_name_xz(p%bb(:,3),idiag_bzmxz)
if (idiag_bx2mxz/=0) call ysum_mn_name_xz(p%bb(:,1)**2,idiag_bx2mxz)
if (idiag_by2mxz/=0) call ysum_mn_name_xz(p%bb(:,2)**2,idiag_by2mxz)
if (idiag_bz2mxz/=0) call ysum_mn_name_xz(p%bb(:,3)**2,idiag_bz2mxz)
if (idiag_bx2mxy/=0) call zsum_mn_name_xy(p%bb(:,1)**2,idiag_bx2mxy)
if (idiag_by2mxy/=0) call zsum_mn_name_xy(p%bb(:,2)**2,idiag_by2mxy)
if (idiag_bz2mxy/=0) call zsum_mn_name_xy(p%bb(:,3)**2,idiag_bz2mxy)
if (idiag_bxbymxy/=0) &
call zsum_mn_name_xy(p%bb(:,1)*p%bb(:,2),idiag_bxbymxy)
if (idiag_bxbzmxy/=0) &
call zsum_mn_name_xy(p%bb(:,1)*p%bb(:,3),idiag_bxbzmxy)
if (idiag_bybzmxy/=0) &
call zsum_mn_name_xy(p%bb(:,2)*p%bb(:,3),idiag_bybzmxy)
if (idiag_bxbymxz/=0) &
call ysum_mn_name_xz(p%bb(:,1)*p%bb(:,2),idiag_bxbymxz)
if (idiag_bxbzmxz/=0) &
call ysum_mn_name_xz(p%bb(:,1)*p%bb(:,3),idiag_bxbzmxz)
if (idiag_bybzmxz/=0) &
call ysum_mn_name_xz(p%bb(:,2)*p%bb(:,3),idiag_bybzmxz)
if (idiag_Exmxz/=0) call ysum_mn_name_xz(p%uxb(:,1),idiag_Exmxz)
if (idiag_Eymxz/=0) call ysum_mn_name_xz(p%uxb(:,2),idiag_Eymxz)
if (idiag_Ezmxz/=0) call ysum_mn_name_xz(p%uxb(:,3),idiag_Ezmxz)
else
!
! idiag_bxmxy and idiag_bymxy also need to be calculated when
! ldiagnos and idiag_bmx and/or idiag_bmy, so
!
! We may need to calculate bxmxy without calculating bmx. The following
! if condition was messing up calculation of bmxy_rms
!
if (ldiagnos) then
if (idiag_bxmxy/=0) call zsum_mn_name_xy(p%bb(:,1),idiag_bxmxy)
if (idiag_bymxy/=0) call zsum_mn_name_xy(p%bb(:,2),idiag_bymxy)
if (idiag_bzmxy/=0) call zsum_mn_name_xy(p%bb(:,3),idiag_bzmxy)
if (idiag_jxmxy/=0) call zsum_mn_name_xy(p%jj(:,1),idiag_jxmxy)
if (idiag_jymxy/=0) call zsum_mn_name_xy(p%jj(:,2),idiag_jymxy)
if (idiag_jzmxy/=0) call zsum_mn_name_xy(p%jj(:,3),idiag_jzmxy)
endif
endif
!
! Debug output.
!
if (headtt .and. lfirst .and. ip<=4) then
call output_pencil(trim(directory)//'/aa.dat',p%aa,3)
call output_pencil(trim(directory)//'/bb.dat',p%bb,3)
call output_pencil(trim(directory)//'/jj.dat',p%jj,3)
call output_pencil(trim(directory)//'/del2A.dat',p%del2a,3)
call output_pencil(trim(directory)//'/JxBr.dat',p%jxbr,3)
call output_pencil(trim(directory)//'/JxB.dat',p%jxb,3)
call output_pencil(trim(directory)//'/df.dat',df(l1:l2,m,n,:),mvar)
endif
!
call timing('daa_dt','finished',mnloop=.true.)
!
endsubroutine daa_dt
!***********************************************************************
subroutine read_magnetic_init_pars(unit,iostat)
!
integer, intent(in) :: unit
integer, intent(inout), optional :: iostat
!
if (present(iostat)) then
read(unit,NML=magnetic_init_pars,ERR=99, IOSTAT=iostat)
else
read(unit,NML=magnetic_init_pars,ERR=99)
endif
!
99 return
!
endsubroutine read_magnetic_init_pars
!***********************************************************************
subroutine write_magnetic_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit,NML=magnetic_init_pars)
!
endsubroutine write_magnetic_init_pars
!***********************************************************************
subroutine read_magnetic_run_pars(unit,iostat)
!
integer, intent(in) :: unit
integer, intent(inout), optional :: iostat
!
if (present(iostat)) then
read(unit,NML=magnetic_run_pars,ERR=99, IOSTAT=iostat)
else
read(unit,NML=magnetic_run_pars,ERR=99)
endif
!
99 return
!
endsubroutine read_magnetic_run_pars
!***********************************************************************
subroutine write_magnetic_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit,NML=magnetic_run_pars)
!
endsubroutine write_magnetic_run_pars
!***********************************************************************
subroutine get_slices_magnetic(f,slices)
!
! Write slices for animation of Magnetic variables.
!
! 26-jul-06/tony: coded
!
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
!
! Loop over slices
!
select case (trim(slices%name))
!
! Magnetic vector potential (code variable)
!
case ('aa')
if (slices%index>=3) then
slices%ready=.false.
else
slices%index=slices%index+1
slices%yz =f(slices%ix,m1:m2 ,n1:n2 ,iax-1+slices%index)
slices%xz =f(l1:l2 ,slices%iy,n1:n2 ,iax-1+slices%index)
slices%xy =f(l1:l2 ,m1:m2 ,slices%iz ,iax-1+slices%index)
slices%xy2=f(l1:l2 ,m1:m2 ,slices%iz2,iax-1+slices%index)
if (slices%index<=3) slices%ready=.true.
endif
!
! Magnetic field (derived variable)
!
case ('bb')
if (slices%index>=3) then
slices%ready=.false.
else
slices%index=slices%index+1
slices%yz =>bb_yz(:,:,slices%index)
slices%xz =>bb_xz(:,:,slices%index)
slices%xy =>bb_xy(:,:,slices%index)
slices%xy2=>bb_xy2(:,:,slices%index)
if (slices%index<=3) slices%ready=.true.
endif
!
! Current density (derived variable)
!
case ('jj')
if (slices%index>=3) then
slices%ready=.false.
else
slices%index=slices%index+1
slices%yz =>jj_yz(:,:,slices%index)
slices%xz =>jj_xz(:,:,slices%index)
slices%xy =>jj_xy(:,:,slices%index)
slices%xy2=>jj_xy2(:,:,slices%index)
if (slices%index<=3) slices%ready=.true.
endif
!
! Magnetic field squared (derived variable)
!
case ('b2')
slices%yz =>b2_yz
slices%xz =>b2_xz
slices%xy =>b2_xy
slices%xy2=>b2_xy2
slices%ready=.true.
!
! Current density (derived variable)
!
case ('jb')
slices%yz =>jb_yz
slices%xz =>jb_xz
slices%xy =>jb_xy
slices%xy2=>jb_xy2
slices%ready=.true.
!
! Plasma beta
!
case ('beta')
slices%yz =>jb_yz
slices%xz =>jb_xz
slices%xy =>jb_xy
slices%xy2=>jb_xy2
slices%ready=.true.
!
endselect
!
endsubroutine get_slices_magnetic
!***********************************************************************
subroutine rprint_magnetic(lreset,lwrite)
!
! Reads and registers print parameters relevant for magnetic fields.
!
! 3-may-02/axel: coded
! 27-may-02/axel: added possibility to reset list
!
use Diagnostics
!
integer :: iname,inamex,inamey,inamez,ixy,ixz,iname_half
logical :: lreset,lwr
logical, optional :: lwrite
!
lwr = .false.
if (present(lwrite)) lwr=lwrite
!
! Reset everything in case of RELOAD.
! (this needs to be consistent with what is defined above!)
!
if (lreset) then
idiag_ab_int=0; idiag_jb_int=0; idiag_b2tm=0; idiag_bjtm=0; idiag_jbtm=0
idiag_b2uzm=0; idiag_b2ruzm=0; idiag_ubbzm=0; idiag_b1m=0; idiag_b2m=0
idiag_bm2=0; idiag_j2m=0; idiag_jm2=0
idiag_abm=0; idiag_abrms=0; idiag_abmh=0
idiag_abumx=0; idiag_abumy=0; idiag_abumz=0
idiag_abmn=0; idiag_abms=0; idiag_jbmh=0; idiag_jbmn=0; idiag_jbms=0
idiag_ajm=0; idiag_cosubm=0; idiag_jbm=0
idiag_uam=0; idiag_ubm=0; idiag_ujm=0
idiag_uxbxm=0; idiag_uybym=0; idiag_uzbzm=0
idiag_fbm=0; idiag_fxbxm=0; idiag_epsM=0; idiag_epsM_LES=0
idiag_epsAD=0; idiag_bxpt=0; idiag_bypt=0; idiag_bzpt=0; idiag_Expt=0
idiag_Eypt=0; idiag_Ezpt=0; idiag_aybym2=0; idiag_exaym2=0
idiag_exjm2=0; idiag_brms=0; idiag_bmax=0; idiag_jrms=0; idiag_jmax=0
idiag_vArms=0; idiag_emag=0; idiag_bxmin=0; idiag_bymin=0; idiag_bzmin=0
idiag_bxmax=0; idiag_bymax=0; idiag_bzmax=0; idiag_vAmax=0; idiag_dtb=0
idiag_arms=0; idiag_amax=0; idiag_Qsm=0; idiag_Qpm=0; idiag_qem=0; idiag_beta1m=0
idiag_beta1max=0; idiag_bxm=0; idiag_bym=0; idiag_bzm=0; idiag_axm=0
idiag_aym=0; idiag_azm=0; idiag_bx2m=0; idiag_by2m=0; idiag_bz2m=0
idiag_bxbymy=0; idiag_bxbzmy=0; idiag_bybzmy=0; idiag_bxbymz=0
idiag_bxbzmz=0; idiag_bybzmz=0; idiag_b2mz=0
idiag_jbmz=0; idiag_abmz=0; idiag_uamz=0
idiag_bxbym=0; idiag_bxbzm=0; idiag_bybzm=0; idiag_djuidjbim=0
idiag_axmz=0; idiag_aymz=0; idiag_azmz=0; idiag_bxmz=0; idiag_bymz=0
idiag_abuxmz=0; idiag_abuymz=0; idiag_abuzmz=0
idiag_uabxmz=0; idiag_uabymz=0; idiag_uabzmz=0
idiag_bzmz=0; idiag_bmx=0; idiag_bmy=0; idiag_bmz=0; idiag_embmz=0
idiag_bmzS2=0; idiag_bmzA2=0
idiag_emxamz3=0; idiag_jmx=0; idiag_jmy=0; idiag_jmz=0; idiag_ambmz=0
idiag_jmbmz=0; idiag_kmz=0; idiag_kx_aa=0
idiag_ambmzh=0;idiag_ambmzn=0;idiag_ambmzs=0; idiag_bmzph=0
idiag_bmzphe=0; idiag_bcosphz=0; idiag_bsinphz=0
idiag_bx2mz=0; idiag_by2mz=0; idiag_bz2mz=0
idiag_bx2rmz=0; idiag_by2rmz=0; idiag_bz2rmz=0
idiag_bxmxy=0; idiag_bymxy=0; idiag_bzmxy=0
idiag_jxmxy=0; idiag_jymxy=0; idiag_jzmxy=0
idiag_bx2mxy=0; idiag_by2mxy=0; idiag_bz2mxy=0; idiag_bxbymxy=0
idiag_bxbzmxy=0; idiag_bybzmxy=0; idiag_bxbymxz=0; idiag_bxbzmxz=0
idiag_bybzmxz=0; idiag_bxmxz=0; idiag_bymxz=0; idiag_bzmxz=0
idiag_axmxz=0; idiag_aymxz=0; idiag_azmxz=0; idiag_Exmxz=0
idiag_Eymxz=0; idiag_Ezmxz=0; idiag_jxbm=0; idiag_uxbm=0; idiag_oxuxbm=0
idiag_jxbxbm=0; idiag_gpxbm=0; idiag_uxDxuxbm=0
idiag_uxbmx=0; idiag_uxbmy=0; idiag_uxbmz=0
idiag_jxbmx=0; idiag_jxbmy=0; idiag_jxbmz=0
idiag_uxbcmx=0; idiag_uxbsmx=0
idiag_uxbcmy=0; idiag_uxbsmy=0; idiag_examz1=0; idiag_examz2=0
idiag_examz3=0; idiag_examx=0; idiag_examy=0; idiag_examz=0
idiag_exjmx=0; idiag_exjmy=0; idiag_exjmz=0; idiag_dexbmx=0
idiag_dexbmy=0; idiag_dexbmz=0; idiag_phibmx=0; idiag_phibmy=0
idiag_phibmz=0; idiag_uxjm=0; idiag_ujxbm=0; idiag_b2divum=0
idiag_b3b21m=0; idiag_b1b32m=0; idiag_b2b13m=0
idiag_udotxbm=0; idiag_uxbdotm=0
idiag_jxbrxm=0; idiag_jxbrym=0; idiag_jxbrzm=0
idiag_jxbr2m=0; idiag_jxbrxmx=0; idiag_jxbrymx=0; idiag_jxbrzmx=0
idiag_jxbrxmy=0; idiag_jxbrymy=0; idiag_jxbrzmy=0; idiag_jxbrxmz=0
idiag_jxbrymz=0; idiag_jxbrzmz=0
idiag_dteta=0; idiag_uxBrms=0; idiag_Bresrms=0
idiag_Rmrms=0; idiag_jfm=0; idiag_va2m=0
idiag_bxmx=0; idiag_bymx=0; idiag_bzmx=0; idiag_bxmy=0
idiag_bymy=0; idiag_bzmy=0; idiag_bx2my=0; idiag_by2my=0; idiag_bz2my=0
idiag_mflux_x=0; idiag_mflux_y=0; idiag_mflux_z=0; idiag_bmxy_rms=0
idiag_brmsh=0; idiag_brmsn=0
idiag_brmss=0; idiag_etatotalmx=0; idiag_etatotalmz=0
idiag_etavamax=0; idiag_etajmax=0; idiag_etaj2max=0; idiag_etajrhomax=0
endif
!
! Check for those quantities that we want to evaluate online.
!
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'ab_int',idiag_ab_int)
call parse_name(iname,cname(iname),cform(iname),'jb_int',idiag_jb_int)
call parse_name(iname,cname(iname),cform(iname),'dteta',idiag_dteta)
call parse_name(iname,cname(iname),cform(iname),'aybym2',idiag_aybym2)
call parse_name(iname,cname(iname),cform(iname),'exaym2',idiag_exaym2)
call parse_name(iname,cname(iname),cform(iname),'exabot',idiag_exabot)
call parse_name(iname,cname(iname),cform(iname),'exatop',idiag_exatop)
call parse_name(iname,cname(iname),cform(iname),'exjm2',idiag_exjm2)
call parse_name(iname,cname(iname),cform(iname),'b2tm',idiag_b2tm)
call parse_name(iname,cname(iname),cform(iname),'bjtm',idiag_bjtm)
call parse_name(iname,cname(iname),cform(iname),'jbtm',idiag_jbtm)
call parse_name(iname,cname(iname),cform(iname),'abm',idiag_abm)
call parse_name(iname,cname(iname),cform(iname),'abmn',idiag_abmn)
call parse_name(iname,cname(iname),cform(iname),'abms',idiag_abms)
call parse_name(iname,cname(iname),cform(iname),'abrms',idiag_abrms)
call parse_name(iname,cname(iname),cform(iname),'abumx',idiag_abumx)
call parse_name(iname,cname(iname),cform(iname),'abumy',idiag_abumy)
call parse_name(iname,cname(iname),cform(iname),'abumz',idiag_abumz)
call parse_name(iname,cname(iname),cform(iname),'ajm',idiag_ajm)
call parse_name(iname,cname(iname),cform(iname),'jbm',idiag_jbm)
call parse_name(iname,cname(iname),cform(iname),'jbmn',idiag_jbmn)
call parse_name(iname,cname(iname),cform(iname),'jbms',idiag_jbms)
call parse_name(iname,cname(iname),cform(iname),'ubm',idiag_ubm)
call parse_name(iname,cname(iname),cform(iname),'uxbxm',idiag_uxbxm)
call parse_name(iname,cname(iname),cform(iname),'uybym',idiag_uybym)
call parse_name(iname,cname(iname),cform(iname),'uzbzm',idiag_uzbzm)
call parse_name(iname,cname(iname),cform(iname),'cosubm',idiag_cosubm)
call parse_name(iname,cname(iname),cform(iname),'uam',idiag_uam)
call parse_name(iname,cname(iname),cform(iname),'ujm',idiag_ujm)
call parse_name(iname,cname(iname),cform(iname),'fbm',idiag_fbm)
call parse_name(iname,cname(iname),cform(iname),'fxbxm',idiag_fxbxm)
call parse_name(iname,cname(iname),cform(iname),'b2ruzm',idiag_b2ruzm)
call parse_name(iname,cname(iname),cform(iname),'b2uzm',idiag_b2uzm)
call parse_name(iname,cname(iname),cform(iname),'ubbzm',idiag_ubbzm)
call parse_name(iname,cname(iname),cform(iname),'b1m',idiag_b1m)
call parse_name(iname,cname(iname),cform(iname),'b2m',idiag_b2m)
call parse_name(iname,cname(iname),cform(iname),'bm2',idiag_bm2)
call parse_name(iname,cname(iname),cform(iname),'j2m',idiag_j2m)
call parse_name(iname,cname(iname),cform(iname),'jm2',idiag_jm2)
call parse_name(iname,cname(iname),cform(iname),'epsM',idiag_epsM)
call parse_name(iname,cname(iname),cform(iname),&
'epsM_LES',idiag_epsM_LES)
call parse_name(iname,cname(iname),cform(iname),'epsAD',idiag_epsAD)
call parse_name(iname,cname(iname),cform(iname),'emag',idiag_emag)
call parse_name(iname,cname(iname),cform(iname),'brms',idiag_brms)
call parse_name(iname,cname(iname),cform(iname),'brmsn',idiag_brmsn)
call parse_name(iname,cname(iname),cform(iname),'brmss',idiag_brmss)
call parse_name(iname,cname(iname),cform(iname),'bmax',idiag_bmax)
call parse_name(iname,cname(iname),cform(iname),'bxmin',idiag_bxmin)
call parse_name(iname,cname(iname),cform(iname),'bymin',idiag_bymin)
call parse_name(iname,cname(iname),cform(iname),'bzmin',idiag_bzmin)
call parse_name(iname,cname(iname),cform(iname),'bxmax',idiag_bxmax)
call parse_name(iname,cname(iname),cform(iname),'bymax',idiag_bymax)
call parse_name(iname,cname(iname),cform(iname),'bzmax',idiag_bzmax)
call parse_name(iname,cname(iname),cform(iname),'jrms',idiag_jrms)
call parse_name(iname,cname(iname),cform(iname),'jmax',idiag_jmax)
call parse_name(iname,cname(iname),cform(iname),'axm',idiag_axm)
call parse_name(iname,cname(iname),cform(iname),'aym',idiag_aym)
call parse_name(iname,cname(iname),cform(iname),'azm',idiag_azm)
call parse_name(iname,cname(iname),cform(iname),'arms',idiag_arms)
call parse_name(iname,cname(iname),cform(iname),'amax',idiag_amax)
call parse_name(iname,cname(iname),cform(iname),'vArms',idiag_vArms)
call parse_name(iname,cname(iname),cform(iname),'vAmax',idiag_vAmax)
call parse_name(iname,cname(iname),cform(iname),'vA2m',idiag_vA2m)
call parse_name(iname,cname(iname),cform(iname),'Qsm',idiag_Qsm)
call parse_name(iname,cname(iname),cform(iname),'Qpm',idiag_Qpm)
call parse_name(iname,cname(iname),cform(iname),'qem',idiag_qem)
call parse_name(iname,cname(iname),cform(iname),&
'beta1m',idiag_beta1m)
call parse_name(iname,cname(iname),cform(iname),&
'beta1max',idiag_beta1max)
call parse_name(iname,cname(iname),cform(iname),'dtb',idiag_dtb)
call parse_name(iname,cname(iname),cform(iname),'bxm',idiag_bxm)
call parse_name(iname,cname(iname),cform(iname),'bym',idiag_bym)
call parse_name(iname,cname(iname),cform(iname),'bzm',idiag_bzm)
call parse_name(iname,cname(iname),cform(iname),'bx2m',idiag_bx2m)
call parse_name(iname,cname(iname),cform(iname),'by2m',idiag_by2m)
call parse_name(iname,cname(iname),cform(iname),'bz2m',idiag_bz2m)
call parse_name(iname,cname(iname),cform(iname),'bxbym',idiag_bxbym)
call parse_name(iname,cname(iname),cform(iname),'bxbzm',idiag_bxbzm)
call parse_name(iname,cname(iname),cform(iname),'bybzm',idiag_bybzm)
call parse_name(iname,cname(iname),cform(iname),&
'djuidjbim',idiag_djuidjbim)
call parse_name(iname,cname(iname),cform(iname),'jxbrxm',idiag_jxbrxm)
call parse_name(iname,cname(iname),cform(iname),'jxbrym',idiag_jxbrym)
call parse_name(iname,cname(iname),cform(iname),'jxbrzm',idiag_jxbrzm)
call parse_name(iname,cname(iname),cform(iname),'jxbr2m',idiag_jxbr2m)
call parse_name(iname,cname(iname),cform(iname),'jxbm',idiag_jxbm)
call parse_name(iname,cname(iname),cform(iname),'uxbm',idiag_uxbm)
call parse_name(iname,cname(iname),cform(iname),'uxbmx',idiag_uxbmx)
call parse_name(iname,cname(iname),cform(iname),'uxbmy',idiag_uxbmy)
call parse_name(iname,cname(iname),cform(iname),'uxbmz',idiag_uxbmz)
call parse_name(iname,cname(iname),cform(iname),'jxbmx',idiag_jxbmx)
call parse_name(iname,cname(iname),cform(iname),'jxbmy',idiag_jxbmy)
call parse_name(iname,cname(iname),cform(iname),'jxbmz',idiag_jxbmz)
call parse_name(iname,cname(iname),cform(iname),'uxbcmx',idiag_uxbcmx)
call parse_name(iname,cname(iname),cform(iname),'uxbcmy',idiag_uxbcmy)
call parse_name(iname,cname(iname),cform(iname),'uxbsmx',idiag_uxbsmx)
call parse_name(iname,cname(iname),cform(iname),'uxbsmy',idiag_uxbsmy)
call parse_name(iname,cname(iname),cform(iname),'examx',idiag_examx)
call parse_name(iname,cname(iname),cform(iname),'examy',idiag_examy)
call parse_name(iname,cname(iname),cform(iname),'examz',idiag_examz)
call parse_name(iname,cname(iname),cform(iname),'exjmx',idiag_exjmx)
call parse_name(iname,cname(iname),cform(iname),'exjmy',idiag_exjmy)
call parse_name(iname,cname(iname),cform(iname),'exjmz',idiag_exjmz)
call parse_name(iname,cname(iname),cform(iname),'dexbmx',idiag_dexbmx)
call parse_name(iname,cname(iname),cform(iname),'dexbmy',idiag_dexbmy)
call parse_name(iname,cname(iname),cform(iname),'dexbmz',idiag_dexbmz)
call parse_name(iname,cname(iname),cform(iname),'phibmx',idiag_phibmx)
call parse_name(iname,cname(iname),cform(iname),'phibmy',idiag_phibmy)
call parse_name(iname,cname(iname),cform(iname),'phibmz',idiag_phibmz)
call parse_name(iname,cname(iname),cform(iname),'uxjm',idiag_uxjm)
call parse_name(iname,cname(iname),cform(iname),'ujxbm',idiag_ujxbm)
call parse_name(iname,cname(iname),cform(iname),'b2divum',idiag_b2divum)
call parse_name(iname,cname(iname),cform(iname),'jxbxbm',idiag_jxbxbm)
call parse_name(iname,cname(iname),cform(iname),'oxuxbm',idiag_oxuxbm)
call parse_name(iname,cname(iname),cform(iname),'gpxbm',idiag_gpxbm)
call parse_name(iname,cname(iname),cform(iname),&
'uxDxuxbm',idiag_uxDxuxbm)
call parse_name(iname,cname(iname),cform(iname),'b3b21m',idiag_b3b21m)
call parse_name(iname,cname(iname),cform(iname),'b1b32m',idiag_b1b32m)
call parse_name(iname,cname(iname),cform(iname),'b2b13m',idiag_b2b13m)
call parse_name(iname,cname(iname),cform(iname),'udotxbm',idiag_udotxbm)
call parse_name(iname,cname(iname),cform(iname),'uxbdotm',idiag_uxbdotm)
call parse_name(iname,cname(iname),cform(iname),'bmx',idiag_bmx)
call parse_name(iname,cname(iname),cform(iname),'bmy',idiag_bmy)
call parse_name(iname,cname(iname),cform(iname),'bmz',idiag_bmz)
call parse_name(iname,cname(iname),cform(iname),'bmzS2',idiag_bmzS2)
call parse_name(iname,cname(iname),cform(iname),'bmzA2',idiag_bmzA2)
call parse_name(iname,cname(iname),cform(iname),'jmx',idiag_jmx)
call parse_name(iname,cname(iname),cform(iname),'jmy',idiag_jmy)
call parse_name(iname,cname(iname),cform(iname),'jmz',idiag_jmz)
call parse_name(iname,cname(iname),cform(iname),'bmzph',idiag_bmzph)
call parse_name(iname,cname(iname),cform(iname),'bmzphe',idiag_bmzphe)
call parse_name(iname,cname(iname),cform(iname),'bcosphz',idiag_bcosphz)
call parse_name(iname,cname(iname),cform(iname),'bsinphz',idiag_bsinphz)
call parse_name(iname,cname(iname),cform(iname),'emxamz3',idiag_emxamz3)
call parse_name(iname,cname(iname),cform(iname),'embmz',idiag_embmz)
call parse_name(iname,cname(iname),cform(iname),'ambmz',idiag_ambmz)
call parse_name(iname,cname(iname),cform(iname),'ambmzn',idiag_ambmzn)
call parse_name(iname,cname(iname),cform(iname),'ambmzs',idiag_ambmzs)
call parse_name(iname,cname(iname),cform(iname),'jmbmz',idiag_jmbmz)
call parse_name(iname,cname(iname),cform(iname),'kmz',idiag_kmz)
call parse_name(iname,cname(iname),cform(iname),'kx_aa',idiag_kx_aa)
call parse_name(iname,cname(iname),cform(iname),'bxpt',idiag_bxpt)
call parse_name(iname,cname(iname),cform(iname),'bypt',idiag_bypt)
call parse_name(iname,cname(iname),cform(iname),'bzpt',idiag_bzpt)
call parse_name(iname,cname(iname),cform(iname),'Expt',idiag_Expt)
call parse_name(iname,cname(iname),cform(iname),'Eypt',idiag_Eypt)
call parse_name(iname,cname(iname),cform(iname),'Ezpt',idiag_Ezpt)
call parse_name(iname,cname(iname),cform(iname),'uxBrms',idiag_uxBrms)
call parse_name(iname,cname(iname),cform(iname),'Bresrms',idiag_Bresrms)
call parse_name(iname,cname(iname),cform(iname),'Rmrms',idiag_Rmrms)
call parse_name(iname,cname(iname),cform(iname),'jfm',idiag_jfm)
call parse_name(iname,cname(iname),cform(iname),'bmxy_rms',idiag_bmxy_rms)
call parse_name(iname,cname(iname),cform(iname),'etasmagm',idiag_etasmagm)
call parse_name(iname,cname(iname),cform(iname),'etasmagmin',idiag_etasmagmin)
call parse_name(iname,cname(iname),cform(iname),'etasmagmax',idiag_etasmagmax)
call parse_name(iname,cname(iname),cform(iname),'etavamax',idiag_etavamax)
call parse_name(iname,cname(iname),cform(iname),'etajmax',idiag_etajmax)
call parse_name(iname,cname(iname),cform(iname),'etaj2max',idiag_etaj2max)
call parse_name(iname,cname(iname),cform(iname),'etajrhomax',idiag_etajrhomax)
call parse_name(iname,cname(iname),cform(iname),'cosjbm',idiag_cosjbm)
call parse_name(iname,cname(iname),cform(iname),'jparallelm',idiag_jparallelm)
call parse_name(iname,cname(iname),cform(iname),'jperpm',idiag_jperpm)
enddo
!
! Quantities which are averaged over half (north-south) the box.
!
iname_half=name_half_max
!
! Magnetic helicity (north and south) of total field.
!
if ((idiag_abmn/=0).or.(idiag_abms/=0)) then
iname_half=iname_half+1
idiag_abmh=iname_half
endif
!
! Current helicity (north and south) of total field.
!
if ((idiag_jbmn/=0).or.(idiag_jbms/=0)) then
iname_half=iname_half+1
idiag_jbmh=iname_half
endif
!
! Magnetic energy (north and south) of total field.
!
if ((idiag_brmsn/=0).or.(idiag_brmss/=0)) then
iname_half=iname_half+1
idiag_brmsh=iname_half
endif
!
! Magnetic helicity (north and south) of mean field.
!
if ((idiag_ambmzn/=0).or.(idiag_ambmzs/=0)) then
iname_half=iname_half+1
idiag_ambmzh=iname_half
endif
!
! Update name_half_max.
!
name_half_max=iname_half
!
! Currently need to force zaverage calculation at every lout step for
! bmx and bmy and bmxy_rms.
!
if ((idiag_bmx+idiag_bmy+idiag_bmxy_rms)>0) ldiagnos_need_zaverages=.true.
!
! Check for those quantities for which we want xy-averages.
!
do inamex=1,nnamex
call parse_name(inamex,cnamex(inamex),cformx(inamex),'bxmx',idiag_bxmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'bymx',idiag_bymx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'bzmx',idiag_bzmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'bx2mx',idiag_bx2mx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'by2mx',idiag_by2mx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'bz2mx',idiag_bz2mx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'bxbymx',idiag_bxbymx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'jxbrxmx',idiag_jxbrxmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'jxbrymx',idiag_jxbrymx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'jxbrzmx',idiag_jxbrzmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'mflux_x',idiag_mflux_x)
call parse_name(inamex,cnamex(inamex),cformx(inamex), &
'etatotalmx',idiag_etatotalmx)
enddo
!
! Check for those quantities for which we want xz-averages.
!
do inamey=1,nnamey
call parse_name(inamey,cnamey(inamey),cformy(inamey),'bxmy',idiag_bxmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey),'bymy',idiag_bymy)
call parse_name(inamey,cnamey(inamey),cformy(inamey),'bzmy',idiag_bzmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'bx2my',idiag_bx2my)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'by2my',idiag_by2my)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'bz2my',idiag_bz2my)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'bxbymy',idiag_bxbymy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'bxbzmy',idiag_bxbzmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'bybzmy',idiag_bybzmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'jxbrxmy',idiag_jxbrxmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'jxbrymy',idiag_jxbrymy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'jxbrzmy',idiag_jxbrzmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey), &
'mflux_y',idiag_mflux_y)
enddo
!
! Check for those quantities for which we want yz-averages.
!
do inamez=1,nnamez
call parse_name(inamez,cnamez(inamez),cformz(inamez),'axmz',idiag_axmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'aymz',idiag_aymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'azmz',idiag_azmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'abuxmz',idiag_abuxmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'abuymz',idiag_abuymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'abuzmz',idiag_abuzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'uabxmz',idiag_uabxmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'uabymz',idiag_uabymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'uabzmz',idiag_uabzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bxmz',idiag_bxmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bymz',idiag_bymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bzmz',idiag_bzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'jxmz',idiag_jxmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'jymz',idiag_jymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'jzmz',idiag_jzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'Exmz',idiag_Exmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'Eymz',idiag_Eymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'Ezmz',idiag_Ezmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bx2mz',idiag_bx2mz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'by2mz',idiag_by2mz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bz2mz',idiag_bz2mz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bx2rmz',idiag_bx2rmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'by2rmz',idiag_by2rmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'bz2rmz',idiag_bz2rmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'bxbymz',idiag_bxbymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'bxbzmz',idiag_bxbzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'bybzmz',idiag_bybzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'b2mz',idiag_b2mz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'jxbrxmz',idiag_jxbrxmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'jxbrymz',idiag_jxbrymz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'jxbrzmz',idiag_jxbrzmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'mflux_z',idiag_mflux_z)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'jbmz',idiag_jbmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'abmz',idiag_abmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'uamz',idiag_uamz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'examz1',idiag_examz1)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'examz2',idiag_examz2)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'examz3',idiag_examz3)
call parse_name(inamez,cnamez(inamez),cformz(inamez), &
'etatotalmz',idiag_etatotalmz)
enddo
!
! Check for those quantities for which we want y-averages.
!
do ixz=1,nnamexz
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'axmxz',idiag_axmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'aymxz',idiag_aymxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'azmxz',idiag_azmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bxmxz',idiag_bxmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bymxz',idiag_bymxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bzmxz',idiag_bzmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bx2mxz',idiag_bx2mxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'by2mxz',idiag_by2mxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bz2mxz',idiag_bz2mxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bxbymxz',idiag_bxbymxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bxbzmxz',idiag_bxbzmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'bybzmxz',idiag_bybzmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'Exmxz',idiag_Exmxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'Eymxz',idiag_Eymxz)
call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'Ezmxz',idiag_Ezmxz)
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)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bymxy',idiag_bymxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bzmxy',idiag_bzmxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'jxmxy',idiag_jxmxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'jymxy',idiag_jymxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'jzmxy',idiag_jzmxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bx2mxy',idiag_bx2mxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'by2mxy',idiag_by2mxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bz2mxy',idiag_bz2mxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bxbymxy',idiag_bxbymxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bxbzmxy',idiag_bxbzmxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'bybzmxy',idiag_bybzmxy)
enddo
!
! Write column, idiag_XYZ, where our variable XYZ is stored.
!
if (lwr) then
write(3,*) 'i_ab_int=',idiag_ab_int
write(3,*) 'i_jb_int=',idiag_jb_int
write(3,*) 'i_dteta=',idiag_dteta
write(3,*) 'i_aybym2=',idiag_aybym2
write(3,*) 'i_exaym2=',idiag_exaym2
write(3,*) 'i_exabot=',idiag_exabot
write(3,*) 'i_exatop=',idiag_exatop
write(3,*) 'i_exjm2=',idiag_exjm2
write(3,*) 'i_b2tm=',idiag_b2tm
write(3,*) 'i_bjtm=',idiag_bjtm
write(3,*) 'i_jbtm=',idiag_jbtm
write(3,*) 'i_abm=',idiag_abm
write(3,*) 'i_abumx=',idiag_abumx
write(3,*) 'i_abumy=',idiag_abumy
write(3,*) 'i_abumz=',idiag_abumz
write(3,*) 'i_abmh=',idiag_abmh
write(3,*) 'i_abmn=',idiag_abmn
write(3,*) 'i_abms=',idiag_abms
write(3,*) 'i_abrms=',idiag_abrms
write(3,*) 'i_ajm=',idiag_ajm
write(3,*) 'i_brmsn=',idiag_brmsn
write(3,*) 'i_brmss=',idiag_brmss
write(3,*) 'i_brmsh=',idiag_brmsh
write(3,*) 'i_jbm=',idiag_jbm
write(3,*) 'i_jbmh=',idiag_abmh
write(3,*) 'i_jbmn=',idiag_abmn
write(3,*) 'i_jbms=',idiag_abms
write(3,*) 'i_ubm=',idiag_ubm
write(3,*) 'i_uxbxm=',idiag_uxbxm
write(3,*) 'i_uybym=',idiag_uybym
write(3,*) 'i_uzbzm=',idiag_uzbzm
write(3,*) 'i_cosubm=',idiag_cosubm
write(3,*) 'i_uam=',idiag_uam
write(3,*) 'i_ujm=',idiag_ujm
write(3,*) 'i_fbm=',idiag_fbm
write(3,*) 'i_fxbxm=',idiag_fxbxm
write(3,*) 'i_b2ruzm=',idiag_b2ruzm
write(3,*) 'i_b2uzm=',idiag_b2uzm
write(3,*) 'i_ubbzm=',idiag_ubbzm
write(3,*) 'i_b1m=',idiag_b1m
write(3,*) 'i_b2m=',idiag_b2m
write(3,*) 'i_bm2=',idiag_bm2
write(3,*) 'i_j2m=',idiag_j2m
write(3,*) 'i_jm2=',idiag_jm2
write(3,*) 'i_epsM=',idiag_epsM
write(3,*) 'i_epsM_LES=',idiag_epsM_LES
write(3,*) 'i_emag=',idiag_emag
write(3,*) 'i_brms=',idiag_brms
write(3,*) 'i_bmax=',idiag_bmax
write(3,*) 'i_jrms=',idiag_jrms
write(3,*) 'i_jmax=',idiag_jmax
write(3,*) 'i_axm=',idiag_axm
write(3,*) 'i_aym=',idiag_aym
write(3,*) 'i_azm=',idiag_azm
write(3,*) 'i_arms=',idiag_arms
write(3,*) 'i_amax=',idiag_amax
write(3,*) 'i_vArms=',idiag_vArms
write(3,*) 'i_vAmax=',idiag_vAmax
write(3,*) 'i_vA2m=',idiag_vA2m
write(3,*) 'i_Qsm=',idiag_Qsm
write(3,*) 'i_Qpm=',idiag_Qpm
write(3,*) 'i_qem=',idiag_qem
write(3,*) 'i_beta1m=',idiag_beta1m
write(3,*) 'i_beta1max=',idiag_beta1max
write(3,*) 'i_dtb=',idiag_dtb
write(3,*) 'i_bxm=',idiag_bxm
write(3,*) 'i_bym=',idiag_bym
write(3,*) 'i_bzm=',idiag_bzm
write(3,*) 'i_bx2m=',idiag_bx2m
write(3,*) 'i_by2m=',idiag_by2m
write(3,*) 'i_bz2m=',idiag_bz2m
write(3,*) 'i_bxbym=',idiag_bxbym
write(3,*) 'i_bxbzm=',idiag_bxbzm
write(3,*) 'i_bybzm=',idiag_bybzm
write(3,*) 'i_djuidjbim=',idiag_djuidjbim
write(3,*) 'i_uxbm=',idiag_uxbm
write(3,*) 'i_jxbm=',idiag_jxbm
write(3,*) 'i_uxbmx=',idiag_uxbmx
write(3,*) 'i_uxbmy=',idiag_uxbmy
write(3,*) 'i_uxbmz=',idiag_uxbmz
write(3,*) 'i_jxbmx=',idiag_jxbmx
write(3,*) 'i_jxbmy=',idiag_jxbmy
write(3,*) 'i_jxbmz=',idiag_jxbmz
write(3,*) 'i_examx=',idiag_examx
write(3,*) 'i_examy=',idiag_examy
write(3,*) 'i_examz=',idiag_examz
write(3,*) 'i_exjmx=',idiag_exjmx
write(3,*) 'i_exjmy=',idiag_exjmy
write(3,*) 'i_exjmz=',idiag_exjmz
write(3,*) 'i_dexbmx=',idiag_dexbmx
write(3,*) 'i_dexbmy=',idiag_dexbmy
write(3,*) 'i_dexbmz=',idiag_dexbmz
write(3,*) 'i_phibmx=',idiag_phibmx
write(3,*) 'i_phibmy=',idiag_phibmy
write(3,*) 'i_phibmz=',idiag_phibmz
write(3,*) 'i_uxjm=',idiag_uxjm
write(3,*) 'i_ujxbm=',idiag_ujxbm
write(3,*) 'i_b2divum=',idiag_b2divum
write(3,*) 'i_oxuxbm=',idiag_oxuxbm
write(3,*) 'i_jxbxbm=',idiag_jxbxbm
write(3,*) 'i_gpxbm=',idiag_gpxbm
write(3,*) 'i_uxDxuxbm=',idiag_uxDxuxbm
write(3,*) 'i_b3b21m=',idiag_b3b21m
write(3,*) 'i_b1b32m=',idiag_b1b32m
write(3,*) 'i_b2b13m=',idiag_b2b13m
write(3,*) 'i_udotxbm=',idiag_udotxbm
write(3,*) 'i_uxbdotm=',idiag_uxbdotm
write(3,*) 'i_axmz=',idiag_axmz
write(3,*) 'i_aymz=',idiag_aymz
write(3,*) 'i_azmz=',idiag_azmz
write(3,*) 'i_abuxmz=',idiag_abuxmz
write(3,*) 'i_abuymz=',idiag_abuymz
write(3,*) 'i_abuzmz=',idiag_abuzmz
write(3,*) 'i_uabxmz=',idiag_uabxmz
write(3,*) 'i_uabymz=',idiag_uabymz
write(3,*) 'i_uabzmz=',idiag_uabzmz
write(3,*) 'i_bxmz=',idiag_bxmz
write(3,*) 'i_bymz=',idiag_bymz
write(3,*) 'i_bzmz=',idiag_bzmz
write(3,*) 'i_jxmz=',idiag_jxmz
write(3,*) 'i_jymz=',idiag_jymz
write(3,*) 'i_jzmz=',idiag_jzmz
write(3,*) 'i_Exmz=',idiag_Exmz
write(3,*) 'i_Eymz=',idiag_Eymz
write(3,*) 'i_Ezmz=',idiag_Ezmz
write(3,*) 'i_bx2mz=',idiag_bx2mz
write(3,*) 'i_by2mz=',idiag_by2mz
write(3,*) 'i_bz2mz=',idiag_bz2mz
write(3,*) 'i_bx2rmz=',idiag_bx2rmz
write(3,*) 'i_by2rmz=',idiag_by2rmz
write(3,*) 'i_bz2rmz=',idiag_bz2rmz
write(3,*) 'i_bxbymz=',idiag_bxbymz
write(3,*) 'i_b2mz=',idiag_b2mz
write(3,*) 'i_bmx=',idiag_bmx
write(3,*) 'i_bmy=',idiag_bmy
write(3,*) 'i_bmz=',idiag_bmz
write(3,*) 'i_bmzS2=',idiag_bmzS2
write(3,*) 'i_bmzA2=',idiag_bmzA2
write(3,*) 'i_jmx=',idiag_jmx
write(3,*) 'i_jmy=',idiag_jmy
write(3,*) 'i_jmz=',idiag_jmz
write(3,*) 'i_bmzph=',idiag_bmzph
write(3,*) 'i_bmzphe=',idiag_bmzphe
write(3,*) 'i_bcosphz=',idiag_bcosphz
write(3,*) 'i_bsinphz=',idiag_bsinphz
write(3,*) 'i_emxamz3=',idiag_emxamz3
write(3,*) 'i_embmz=',idiag_embmz
write(3,*) 'i_ambmz=',idiag_ambmz
write(3,*) 'i_ambmzh=',idiag_ambmzh
write(3,*) 'i_ambmzn=',idiag_ambmzn
write(3,*) 'i_ambmzs=',idiag_ambmzs
write(3,*) 'i_jmbmz=',idiag_jmbmz
write(3,*) 'i_kmz=',idiag_kmz
write(3,*) 'i_kx_aa=',idiag_kx_aa
write(3,*) 'i_bxpt=',idiag_bxpt
write(3,*) 'i_bypt=',idiag_bypt
write(3,*) 'i_bzpt=',idiag_bzpt
write(3,*) 'i_Expt=',idiag_Expt
write(3,*) 'i_Eypt=',idiag_Eypt
write(3,*) 'i_Ezpt=',idiag_Ezpt
write(3,*) 'i_bxmxy=',idiag_bxmxy
write(3,*) 'i_bymxy=',idiag_bymxy
write(3,*) 'i_bzmxy=',idiag_bzmxy
write(3,*) 'i_jxmxy=',idiag_jxmxy
write(3,*) 'i_jymxy=',idiag_jymxy
write(3,*) 'i_jzmxy=',idiag_jzmxy
write(3,*) 'i_axmxz=',idiag_axmxz
write(3,*) 'i_aymxz=',idiag_aymxz
write(3,*) 'i_azmxz=',idiag_azmxz
write(3,*) 'i_bxmxz=',idiag_bxmxz
write(3,*) 'i_bymxz=',idiag_bymxz
write(3,*) 'i_bzmxz=',idiag_bzmxz
write(3,*) 'i_bx2mxz=',idiag_bx2mxz
write(3,*) 'i_by2mxz=',idiag_by2mxz
write(3,*) 'i_bz2mxz=',idiag_bz2mxz
write(3,*) 'i_Exmxz=',idiag_Exmxz
write(3,*) 'i_Eymxz=',idiag_Eymxz
write(3,*) 'i_Ezmxz=',idiag_Ezmxz
write(3,*) 'i_jbmz=',idiag_jbmz
write(3,*) 'i_abmz=',idiag_abmz
write(3,*) 'i_uamz=',idiag_uamz
write(3,*) 'i_examz1=',idiag_examz1
write(3,*) 'i_examz2=',idiag_examz2
write(3,*) 'i_examz3=',idiag_examz3
write(3,*) 'i_bx2mxy=',idiag_bx2mxy
write(3,*) 'i_by2mxy=',idiag_by2mxy
write(3,*) 'i_bz2mxy=',idiag_bz2mxy
write(3,*) 'i_bxbymxy=',idiag_bxbymxy
write(3,*) 'i_bxbzmxy=',idiag_bxbzmxy
write(3,*) 'i_bybzmxy=',idiag_bybzmxy
write(3,*) 'i_bxbymxz=',idiag_bxbymxz
write(3,*) 'i_bxbzmxz=',idiag_bxbzmxz
write(3,*) 'i_bybzmxz=',idiag_bybzmxz
write(3,*) 'i_uxBrms=',idiag_uxBrms
write(3,*) 'i_Bresrms=',idiag_Bresrms
write(3,*) 'i_Rmrms=',idiag_Rmrms
write(3,*) 'i_jfm=',idiag_jfm
write(3,*) 'i_bxmx=',idiag_bxmx
write(3,*) 'i_bxmy=',idiag_bxmy
write(3,*) 'i_bymy=',idiag_bymy
write(3,*) 'i_bzmy=',idiag_bzmy
write(3,*) 'i_bx2my=',idiag_bx2my
write(3,*) 'i_by2my=',idiag_by2my
write(3,*) 'i_bz2my=',idiag_bz2my
write(3,*) 'i_mflux_x=',idiag_mflux_x
write(3,*) 'i_mflux_y=',idiag_mflux_y
write(3,*) 'i_mflux_z=',idiag_mflux_z
write(3,*) 'nname=',nname
write(3,*) 'nnamexy=',nnamexy
write(3,*) 'nnamexz=',nnamexz
write(3,*) 'nnamez=',nnamez
write(3,*) 'iaa=',iaa
write(3,*) 'iax=',iax
write(3,*) 'iay=',iay
write(3,*) 'iaz=',iaz
write(3,*) 'idiag_bmxy_rms=',idiag_bmxy_rms
write(3,*) 'idiag_cosjbm=',idiag_cosjbm
write(3,*) 'i_etavamax=', idiag_etavamax
write(3,*) 'i_etajmax=', idiag_etajmax
write(3,*) 'i_etaj2max=', idiag_etaj2max
write(3,*) 'i_etajrhomax=', idiag_etajrhomax
write(3,*) 'idiag_jparallel=',idiag_jparallelm
write(3,*) 'idiag_jperp=',idiag_jperpm
endif
!
endsubroutine rprint_magnetic
!***********************************************************************
endmodule Magnetic