! $Id$ ! ! 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 2 ! MAUX CONTRIBUTION 0 ! ! PENCILS PROVIDED aa(3); a2; aij(3,3); bb(3); bbb(3); ab; uxb(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); beta1; uga(3); djuidjbi; jo ! PENCILS PROVIDED ujxb; oxu(3); oxuxb(3); jxbxb(3); jxbrxb(3) ! PENCILS PROVIDED glnrhoxb(3); del4a(3); del6a(3); oxj(3); diva ! PENCILS PROVIDED jij(3,3); sj; ss12; mf_EMF(3); mf_EMFdotB ! PENCILS PROVIDED cosjb,jparallel;jperp ! PENCILS PROVIDED cosub ! !*************************************************************** module Magnetic ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! 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,ny,3) :: bb_xy3,jj_xy3 real, target, dimension (nx,ny,3) :: bb_xy4,jj_xy4 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 (nx,ny) :: b2_xy3,jb_xy3 real, target, dimension (nx,ny) :: b2_xy4,jb_xy4 real, target, dimension (ny,nz) :: b2_yz,jb_yz real, target, dimension (nx,nz) :: b2_xz,jb_xz ! real, target, dimension (nx,ny) :: beta1_xy real, target, dimension (nx,ny) :: beta1_xy2 real, target, dimension (nx,ny) :: beta1_xy3 real, target, dimension (nx,ny) :: beta1_xy4 real, target, dimension (ny,nz) :: beta1_yz real, target, dimension (nx,nz) :: beta1_xz ! real, dimension (mx,my) :: alpha_input ! ! Parameters ! integer, parameter :: nresi_max=4 ! real, dimension (ninit) :: amplaa=0.0,kx_aa=1.,ky_aa=1.,kz_aa=1. real, dimension (ninit) :: phasex_aa=0., phasey_aa=0., phasez_aa=0. character (len=labellen), dimension(ninit) :: initaa='nothing' character (len=labellen) :: borderaa='nothing' character (len=labellen), dimension(nresi_max) :: iresistivity='' character (len=labellen) :: alpha_xprofile='const',alpha_yprofile='const' character (len=labellen) :: Omega_xprofile='const',Omega_yprofile='const' character (len=labellen) :: fring_profile='tanh' ! input parameters complex, dimension(3) :: coefaa=(/0.,0.,0./), coefbb=(/0.,0.,0./) real, dimension(3) :: B_ext=(/0.,0.,0./),B1_ext,B_ext_tmp,eta_aniso_hyper3 real, dimension(3) :: axisr1=(/0,0,1/),dispr1=(/0.,0.5,0./) real, dimension(3) :: axisr2=(/1,0,0/),dispr2=(/0.,-0.5,0./) real, dimension(3) :: axisr3=(/1,0,0/),dispr3=(/0.,-0.5,0./) ! ! profile functions ! real, dimension(nx) :: alpha_x,dalpha_x,Omega_x,dOmega_x real, dimension(my) :: alpha_y,dalpha_y,Omega_y,dOmega_y,Omega_tmp0 ! real, target :: zmode=1. !(temporary) ! ! profile parameters ! real :: alpha_x1=0.,alpha_dx1=0. real :: Omega_x1=0.,Omega_dx1=0. real :: Omega_yc1=0.,Omega_yc2=0. ! real :: fring1=0.,Iring1=0.,Rring1=1.,wr1=0.3 real :: fring2=0.,Iring2=0.,Rring2=1.,wr2=0.3 real :: fring3=0.,Iring3=0.,Rring3=1.,wr3=0.3 real :: radius=.1,epsilonaa=1e-2,widthaa=.5,x0aa=0.,z0aa=0. real :: by_left=0.,by_right=0.,bz_left=0.,bz_right=0. real :: ABC_A=1.,ABC_B=1.,ABC_C=1. real :: bthresh=0.,bthresh_per_brms=0.,brms=0.,bthresh_scl=1. real :: eta_shock=0. real :: rhomin_jxb=0.,va2max_jxb=0. real :: omega_Bz_ext=0. real :: mu_r=-0.5 !(still needed for backwards compatibility) real :: mu_ext_pot=-0.5,inclaa=0. real :: mu012=.5 !(=1/2mu0) real :: rescale_aa=0. real :: ampl_B0=0.,D_smag=0.17,B_ext21,B_ext11 real :: Omega_ampl=0. real :: rmode=1.,rm_int=0.,rm_ext=0. real :: alpha_effect=0.,alpha_quenching=0.,delta_effect=0.,meanfield_etat=0. real :: alpha_eps=0. real :: alpha_equator=impossible,alpha_equator_gap=0.,alpha_gap_step=0. real :: alpha_cutoff_up=0.,alpha_cutoff_down=0. real :: meanfield_Qs=1.,meanfield_Qp=1. real :: meanfield_Bs=1.,meanfield_Bp=1. real :: meanfield_kf=1.,meanfield_etaB=0. real :: displacement_gun=0. real :: pertamplaa=0. real :: initpower_aa=0.,cutoff_aa=0.,brms_target=1.,rescaling_fraction=1. real :: phase_beltrami=0., ampl_beltrami=0. real :: bmz=0, bmz_beltrami_phase=0. real :: taareset=0.,daareset=0. real :: center1_x=0., center1_y=0., center1_z=0. real :: fluxtube_border_width=impossible 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 :: llorentzforce=.true.,linduction=.true. logical :: lalpha_phi_equation=.true. logical :: lresi_eta_const=.false. logical :: lresi_etaSS=.false. logical :: lresi_hyper2=.false. logical :: lresi_hyper3=.false. logical :: lresi_hyper3_polar=.false. logical :: lresi_hyper3_strict=.false. logical :: lresi_zdep=.false., lresi_dust=.false., lresi_rdep=.false. logical :: lresi_hyper3_aniso=.false. logical :: lresi_eta_shock=.false. logical :: lresi_eta_shock_perp=.false. logical :: lresi_shell=.false. logical :: lresi_smagorinsky=.false. logical :: lresi_smagorinsky_cross=.false. logical, target, dimension (3) :: lfrozen_bb_bot=(/.false.,.false.,.false./) logical, target, dimension (3) :: lfrozen_bb_top=(/.false.,.false.,.false./) logical :: reinitialize_aa=.false., lohmic_heat=.true., lneutralion_heat=.true. logical :: lB_ext_pot=.false. logical :: lforce_free_test=.false. logical :: lmeanfield_theory=.false.,lOmega_effect=.false. logical :: lmeanfield_noalpm=.false. logical :: lmeanfield_jxb=.false.,lmeanfield_jxb_with_vA2=.false. logical :: lgauss=.false. logical :: lbb_as_aux=.false.,ljj_as_aux=.false. logical :: lbbt_as_aux=.false.,ljjt_as_aux=.false. logical :: lbext_curvilinear=.true., lcheck_positive_va2=.false. logical :: lreset_aa=.false. character (len=labellen) :: pertaa='zero' ! 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,lbext_curvilinear, & lbbt_as_aux,ljjt_as_aux, lneutralion_heat, & center1_x,center1_y,center1_z,fluxtube_border_width, & va2max_jxb,va2power_jxb ! ! run parameters real :: eta=0.,eta1=0.,eta_hyper2=0.,eta_hyper3=0.,height_eta=0.,eta_out=0. real :: meanfield_molecular_eta=0. real :: eta_int=0.,eta_ext=0.,wresistivity=.01 real :: tau_aa_exterior=0. real :: sigma_ratio=1.,eta_width=0.,eta_z0=1. real :: alphaSSm=0. real :: k1_ff=1.,ampl_ff=1.,swirl=1. real :: k1x_ff=1.,k1y_ff=1.,k1z_ff=1. real :: inertial_length=0.,linertial_2 real :: forcing_continuous_aa_phasefact=1. real :: forcing_continuous_aa_amplfact=1., ampl_fcont_aa=1. real :: LLambda_aa=0. real, dimension(nx,my) :: eta_r real, dimension(nx,my,3) :: geta_r real, dimension(mz) :: coskz,sinkz,eta_z real, dimension(mz,3) :: geta_z logical :: lfreeze_aint=.false.,lfreeze_aext=.false. logical :: lweyl_gauge=.false. logical :: lupw_aa=.false. logical :: lforcing_cont_aa=.false. logical :: lelectron_inertia=.false. logical :: lkinematic=.false. logical :: luse_Bext_in_b2=.false. logical :: lmean_friction=.false. character (len=labellen) :: zdep_profile='fs' character (len=labellen) :: rdep_profile='schnack89' character (len=labellen) :: iforcing_continuous_aa='fixed_swirl' ! namelist /magnetic_run_pars/ & eta,eta1,eta_hyper2,eta_hyper3,B_ext,omega_Bz_ext, & lmeanfield_theory,alpha_effect,alpha_quenching,delta_effect, & alpha_eps, & alpha_xprofile,alpha_x1,alpha_dx1, & alpha_yprofile, & lalpha_phi_equation, & lOmega_effect,Omega_ampl, & Omega_xprofile,Omega_x1,Omega_dx1, & Omega_yprofile,Omega_yc1,Omega_yc2, & lmeanfield_noalpm, & meanfield_etat, lohmic_heat, & lmeanfield_jxb,lmeanfield_jxb_with_vA2, & meanfield_Qs, meanfield_Qp, & meanfield_Bs, meanfield_Bp, & meanfield_kf,meanfield_etaB, & alpha_equator,alpha_equator_gap,alpha_gap_step,& alpha_cutoff_up,alpha_cutoff_down,& height_eta,eta_out,tau_aa_exterior, & kx_aa,ky_aa,kz_aa,ABC_A,ABC_B,ABC_C, & lforcing_cont_aa,iforcing_continuous_aa, & forcing_continuous_aa_phasefact, & forcing_continuous_aa_amplfact, & k1_ff,ampl_ff,swirl,radius, & k1x_ff,k1y_ff,k1z_ff,lcheck_positive_va2, & lmean_friction,LLambda_aa, & bthresh,bthresh_per_brms, & iresistivity,lweyl_gauge,lupw_aa, & alphaSSm, & eta_int,eta_ext,eta_shock,wresistivity, & rhomin_jxb,va2max_jxb,va2power_jxb,llorentzforce,linduction, & reinitialize_aa,rescale_aa,lB_ext_pot, & displacement_gun, & pertaa,pertamplaa,D_smag,brms_target,rescaling_fraction, & lfreeze_aint,lfreeze_aext, & sigma_ratio,zdep_profile,eta_width,eta_z0, & borderaa,eta_aniso_hyper3, & lbext_curvilinear, & lbb_as_aux,ljj_as_aux,lkinematic, & lbbt_as_aux,ljjt_as_aux, & lneutralion_heat, lreset_aa, daareset, & luse_Bext_in_b2, ampl_fcont_aa ! ! diagnostic variables (need to be consistent with reset list below) integer :: idiag_aphi2m=0 ! DIAG_DOC: $\left$ integer :: idiag_bphi2m=0 ! DIAG_DOC: $\left$ integer :: idiag_bpol2m=0 ! DIAG_DOC: $\left$ ! ! plus all the others currently in nomagnetic ! integer :: idiag_b2m=0,idiag_bm2=0,idiag_j2m=0,idiag_jm2=0,idiag_abm=0 integer :: idiag_jbm=0,idiag_epsM=0,idiag_vArms=0,idiag_vAmax=0 integer :: idiag_brms=0,idiag_bmax=0,idiag_jrms=0,idiag_jmax=0 integer :: idiag_bx2m=0, idiag_by2m=0, idiag_bz2m=0,idiag_bmz=0 integer :: idiag_bxmz=0,idiag_bymz=0,idiag_bzmz=0,idiag_bmx=0,idiag_bmy=0 integer :: idiag_bxmxy=0,idiag_bymxy=0,idiag_bzmxy=0 integer :: idiag_uxbm=0,idiag_oxuxbm=0,idiag_jxbxbm=0,idiag_uxDxuxbm=0 integer :: idiag_b2mphi=0 integer :: idiag_bmxy_rms=0 integer :: idiag_bsinphz=0 integer :: idiag_bcosphz=0 ! contains ! !*********************************************************************** subroutine register_magnetic() ! ! Initialise variables which should know that we solve for the vector ! potential: iaphi and bphi, etc; increase nvar accordingly ! ! 1-may-02/wolf: coded ! use FArrayManager ! call farray_register_pde('aphi',iaphi) call farray_register_pde('bphi',ibphi) ! ! Identify version number. ! if (lroot) call svn_id( & "$Id$") ! ! Writing files for use with IDL ! if (lroot) then if (maux == 0) then if (nvar < mvar) write(4,*) ',aphi,bphi $' if (nvar == mvar) write(4,*) ',aphi,bphi' else write(4,*) ',aphi,bphi $' endif write(15,*) 'aphi = fltarr(mx,my,mz)*one' write(15,*) 'bphi = fltarr(mx,my,mz)*one' endif ! endsubroutine register_magnetic !*********************************************************************** subroutine initialize_magnetic(f) ! ! 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 BorderProfiles, only: request_border_driving use FArrayManager use SharedVariables use Sub, only: erfunc ! real, dimension (mx,my,mz,mfarray) :: f integer :: i,ierr ! ! 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.) then eta=1./eta1 endif ! ! calculate B_ext21 = 1/B_ext**2 and the unit vector B1_ext = B_ext/|B_ext| ! B_ext21=B_ext(1)**2+B_ext(2)**2 if (B_ext21/=0.) then B_ext21=1./B_ext21 else B_ext21=1. endif B_ext11=sqrt(B_ext21) B1_ext=B_ext*B_ext11 ! ! rescale magnetic field by a factor reinitialize_aa ! if (reinitialize_aa) then f(:,:,:,iax:iaz)=rescale_aa*f(:,:,:,iax:iaz) endif ! ! set lrescaling_magnetic=T if linit_aa=T ! if (lreset_aa) then lrescaling_magnetic=.true. endif ! if (lfreeze_aint) lfreeze_varint(iax:iaz) = .true. if (lfreeze_aext) lfreeze_varext(iax:iaz) = .true. ! ! Initialize resistivity. ! if (iresistivity(1)=='') iresistivity(1)='eta-const' ! default lresi_eta_const=.false. lresi_hyper2=.false. lresi_hyper3=.false. lresi_hyper3_polar=.false. lresi_hyper3_strict=.false. lresi_hyper3_aniso=.false. lresi_eta_shock=.false. lresi_eta_shock_perp=.false. lresi_smagorinsky=.false. lresi_smagorinsky_cross=.false. ! do i=1,nresi_max select case (iresistivity(i)) case ('eta-const') if (lroot) print*, 'resistivity: constant eta' lresi_eta_const=.true. case ('etaSS') if (lroot) print*, 'resistivity: etaSS (Shakura-Sunyaev)' lresi_etaSS=.true. case ('hyper2') if (lroot) print*, 'resistivity: hyper2' lresi_hyper2=.true. case ('hyper3') if (lroot) print*, 'resistivity: hyper3' lresi_hyper3=.true. case ('hyper3_cyl','hyper3-cyl','hyper3_sph','hyper3-sph') if (lroot) print*, 'resistivity: hyper3 curvilinear' lresi_hyper3_polar=.true. case ('hyper3_strict') if (lroot) print*, 'resistivity: strict hyper3 with positive definite heating rate' lresi_hyper3_strict=.true. case ('rdep') if (lroot) print*, 'resistivity: r-dependent' lresi_rdep=.true. call eta_rdep(eta_r,geta_r,rdep_profile) case ('zdep') if (lroot) print*, 'resistivity: z-dependent' lresi_zdep=.true. call eta_zdep(eta_z,geta_z,zdep_profile) case ('dust') if (lroot) print*, 'resistivity: depending on dust density' lresi_dust=.true. case ('hyper3-aniso') if (lroot) print*, 'resistivity: hyper3_aniso' lresi_hyper3_aniso=.true. case ('shell') if (lroot) print*, 'resistivity: shell' lresi_shell=.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 ('shock-perp') if (lroot) print*, 'resistivity: shock_perp' lresi_eta_shock_perp=.true. if (.not. lshock) & call fatal_error('initialize_magnetic', & 'shock resistivity, but module setting SHOCK=noshock') case ('smagorinsky') if (lroot) print*, 'resistivity: smagorinsky' lresi_smagorinsky=.true. case ('smagorinsky-cross') if (lroot) print*, 'resistivity: smagorinsky_cross' lresi_smagorinsky_cross=.true. 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.and.meanfield_etat==0.0)) & call warning('initialize_magnetic', & 'Resistivity coefficient eta is zero!') if (lresi_hyper2.and.eta_hyper2==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_hyper2 is zero!') if (lresi_hyper3.and.eta_hyper3==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_hyper3 is zero!') if (lresi_hyper3_polar.and.eta_hyper3==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_hyper3 is zero!') if (lresi_hyper3_strict.and.eta_hyper3==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_hyper3 is zero!') if ( (lresi_hyper3_aniso) .and. & ((eta_aniso_hyper3(1)==0. .and. nxgrid/=1 ).or. & (eta_aniso_hyper3(2)==0. .and. nygrid/=1 ).or. & (eta_aniso_hyper3(3)==0. .and. nzgrid/=1 )) ) & call fatal_error('initialize_magnetic', & 'A resistivity coefficient of eta_aniso_hyper3 is zero!') if (lresi_eta_shock.and.eta_shock==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_shock is zero!') if (lresi_eta_shock_perp.and.eta_shock==0.0) & call fatal_error('initialize_magnetic', & 'Resistivity coefficient eta_shock is zero!') ! if (lentropy .and. lohmic_heat .and. .not. lresi_eta_const) & ! call fatal_error('initialize_magnetic', & ! 'Resistivity heating only works with regular resistivity!') if (lresi_hyper2.and.lresi_hyper3) & call warning('initialize_magnetic', & '4th & 6th order hyperdiffusion are both set. ' // & 'Timestep is currently only sensitive to fourth order.') endif ! ! alpha profile ! ============= ! select case (alpha_xprofile) case ('const'); alpha_x=1.; dalpha_x=0. case ('erfx'); alpha_x=.5*(1.+erfunc((x(l1:l2)-alpha_x1)/alpha_dx1)) dalpha_x=exp(-((x(l1:l2)-alpha_x1)/alpha_dx1)**2)/(alpha_dx1*sqrtpi) case default if (lroot) print*,'No such such value for alpha_yprofile' call fatal_error('initialize_magnetic','') endselect ! ! y direction ! select case (alpha_yprofile) case ('const'); alpha_y=1.; dalpha_y=0. case ('cosy'); alpha_y=cos(y); dalpha_y=-sin(y) case ('cosy*sin2y'); alpha_y=1.5*sqrt(3.)*cos(y)*sin(y)**2 dalpha_y=1.5*sqrt(3.)*(2.-3.*sin(y)**2)*sin(y) case ('read') print*,'read alpha profile' open(1,file='alpha_input.dat',form='unformatted') read(1) alpha_y !-- read(1) alpha_input close(1) case default if (lroot) print*,'No such such value for alpha_yprofile' call fatal_error('initialize_magnetic','') endselect ! ! Omega effect and its gradient ! ============================= ! if (lOmega_effect) then select case (Omega_xprofile) case ('const'); Omega_x=1.; dOmega_x=0. case ('Omega=x'); Omega_x=x(l1:l2); dOmega_x=1. case ('erfx'); Omega_x=.5*(1.+erfunc((x(l1:l2)-Omega_x1)/Omega_dx1)) dOmega_x=exp(-((x(l1:l2)-Omega_x1)/Omega_dx1)**2)/(Omega_dx1*sqrtpi) case default if (lroot) print*,'No such such value for Omega_xprofile' call fatal_error('initialize_magnetic','') endselect ! ! y direction; compute Omega_y and dOmega_y = (dOmega_y/dy) ! select case (Omega_yprofile) case ('const'); Omega_y=1.; dOmega_y=0. case ('c1+c2*cos2y'); Omega_y=Omega_yc1+Omega_yc2*cos(y)**2 dOmega_y=-2.*Omega_yc2*cos(y)*sin(y) case default if (lroot) print*,'No such such value for Omega_yprofile' call fatal_error('initialize_magnetic','') endselect endif ! ! if meanfield theory is invoked, we want to send meanfield_etat to ! other subroutines ! call put_shared_variable('lmeanfield_theory',lmeanfield_theory,ierr) if (lmeanfield_theory) then call put_shared_variable('meanfield_etat',meanfield_etat,ierr) call put_shared_variable('eta',eta,ierr) endif ! ! Tell the BorderProfiles module if we intend to use border driving, so ! that the module can request the right pencils. ! if (borderaa/='nothing') call request_border_driving(borderaa) ! ! 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 else if (lroot) print*, 'initialize_magnetic: ibb = ', ibb call farray_index_append('ibb',ibb) call farray_index_append('ibx',ibx) call farray_index_append('iby',iby) call farray_index_append('ibz',ibz) 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 else if (lroot) print*, 'initialize_magnetic: ijj = ', ijj call farray_index_append('ijj',ijj) call farray_index_append('ijx',ijx) call farray_index_append('ijy',ijy) call farray_index_append('ijz',ijz) endif endif ! ! After a reload, we need to rewrite index.pro, but the auxiliary ! arrays are already allocated and must not be allocated again. ! if (lbbt_as_aux) then if (ibbt==0) then call farray_register_auxiliary('bbt',ibbt,vector=3) ibxt=ibbt ibyt=ibbt+1 ibzt=ibbt+2 else if (lroot) print*, 'initialize_velocity: ibbt = ', ibbt call farray_index_append('ibbt',ibbt) call farray_index_append('ibxt',ibxt) call farray_index_append('ibyt',ibyt) call farray_index_append('ibzt',ibzt) endif endif ! if (ljjt_as_aux) then if (ijjt==0) then call farray_register_auxiliary('jjt',ijjt,vector=3) ijxt=ijjt ijyt=ijjt+1 ijzt=ijjt+2 else if (lroot) print*, 'initialize_velocity: ijjt = ', ijjt call farray_index_append('ijjt',ijjt) call farray_index_append('ijxt',ijxt) call farray_index_append('ijyt',ijyt) call farray_index_append('ijzt',ijzt) endif endif ! if (any(initaa=='Alfven-zconst')) then call put_shared_variable('zmode',zmode,ierr) if (ierr/=0) call fatal_error('initialize_magnetic',& 'there was a problem when sharing zmode') endif ! call put_shared_variable('lfrozen_bb_bot',lfrozen_bb_bot,ierr) if (ierr/=0) call fatal_error('initialize_magnetic',& 'there was a problem when sharing lfrozen_bb_bot') call put_shared_variable('lfrozen_bb_top',lfrozen_bb_top,ierr) if (ierr/=0) call fatal_error('initialize_magnetic',& 'there was a problem when sharing lfrozen_bb_top') ! 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 InitialCondition, only: initial_condition_aa use Mpicomm use SharedVariables use Sub ! real, dimension (mx,my,mz,mfarray) :: f ! real, dimension (mz) :: tmp real, dimension (nx,3) :: bb real, dimension (nx) :: b2,fact real :: beq2 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(:,:,:,iaphi) = 0.; f(:,:,:,ibphi) = 0. case ('gaussian-noise'); call gaunoise(amplaa(j),f,iaphi); call gaunoise(amplaa(j),f,ibphi) case ('dipolar-field'); f(:,:,:,ibphi) = 0.; call phi_siny_over_r2(f,iaphi) ! 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) ! ! allow for pressure equilibrium (for isothermal tube) ! assume that ghost zones have already been set. ! corrected expression below for gamma /= 1 case. ! The beq2 expression for 2*mu0*p is not general yet. ! if (lpress_equil.or.lpress_equil_via_ss) then if (lroot) print*,'init_aa: adjust lnrho to have pressure equilib; cs0=',cs0 do n=n1,n2 do m=m1,m2 call curl(f,iaa,bb) call dot2_mn(bb,b2) if (gamma==1.) then f(l1:l2,m,n,ilnrho)=f(l1:l2,m,n,ilnrho)-b2/(2.*cs0**2) else beq2=2.*rho0*cs0**2 fact=max(1e-6,1.-b2/beq2) if (lentropy.and.lpress_equil_via_ss) then f(l1:l2,m,n,iss)=f(l1:l2,m,n,iss)+fact/gamma else f(l1:l2,m,n,ilnrho)=f(l1:l2,m,n,ilnrho)+fact/gamma_m1 endif endif enddo enddo endif ! 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. ! ! jj pencil always needed when in Weyl gauge ! if ( height_eta/=0. .or. ip<=4 .or. & (lweyl_gauge) .or. (lspherical_coords) ) & ! WL: but doesn't seem to be needed for the cylindrical case lpenc_requested(i_jj)=.true. if ((eta/=0..or.meanfield_etat/=0.).and. & (.not.lweyl_gauge)) lpenc_requested(i_del2a)=.true. if (dvid/=0.) lpenc_video(i_jb)=.true. if (lresi_eta_const .or. lresi_shell .or. & lresi_eta_shock .or. lresi_smagorinsky .or. & lresi_zdep .or. lresi_rdep .or. & lresi_smagorinsky_cross) lpenc_requested(i_del2a)=.true. if (lresi_eta_shock) then lpenc_requested(i_shock)=.true. if (lweyl_gauge) then lpenc_requested(i_jj)=.true. else lpenc_requested(i_gshock)=.true. lpenc_requested(i_diva)=.true. endif endif if (lresi_shell) then lpenc_requested(i_r_mn)=.true. lpenc_requested(i_evr)=.true. endif if (lresi_eta_shock_perp) then lpenc_requested(i_shock_perp)=.true. if (lweyl_gauge) then lpenc_requested(i_jj)=.true. else lpenc_requested(i_gshock_perp)=.true. lpenc_requested(i_diva)=.true. endif endif if (lupw_aa) then lpenc_requested(i_uu)=.true. lpenc_requested(i_aij)=.true. endif if (lresi_shell .or. lresi_zdep .or. lresi_rdep) lpenc_requested(i_diva)=.true. if (lresi_smagorinsky_cross) lpenc_requested(i_jo)=.true. if (lresi_hyper2) lpenc_requested(i_del4a)=.true. if (lresi_hyper3) lpenc_requested(i_del6a)=.true. !WL: for the cylindrical case, lpencil_check says graddiva is not needed if (lspherical_coords) lpenc_requested(i_graddiva)=.true. if (lentropy .or. lresi_smagorinsky .or. ltemperature) then lpenc_requested(i_j2)=.true. endif if (lresi_dust) lpenc_requested(i_rhop)=.true. if (lentropy .or. ltemperature .or. ldt) lpenc_requested(i_rho1)=.true. if (lentropy .or. ltemperature) lpenc_requested(i_TT1)=.true. if (ltemperature) lpenc_requested(i_cv1)=.true. ! ! for mean-field modelling ! if (lmeanfield_jxb) then lpenc_requested(i_va2)=.true. lpenc_requested(i_rho1)=.true. endif ! ! ambipolar diffusion ! if (lhydro .and. llorentzforce) & lpenc_requested(i_jxbr)=.true. if (lresi_smagorinsky_cross .or. delta_effect/=0.) & lpenc_requested(i_oo)=.true. if (lmeanfield_theory) then if (meanfield_etat/=0. .or. alpha_effect/=0. .or. delta_effect/=0.) & lpenc_requested(i_mf_EMF)=.true. if (delta_effect/=0.) lpenc_requested(i_oxj)=.true. 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_va2)) then lpencil_in(i_b2)=.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_beta1)) 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 if (.not.lcartesian_coords) lpencil_in(i_aa)=.true. 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_mf_EMFdotB)) then lpencil_in(i_mf_EMF)=.true. lpencil_in(i_bb)=.true. endif if (lpencil_in(i_mf_EMF)) then if (lspherical_coords) then lpencil_in(i_jj)=.true. lpencil_in(i_graddivA)=.true. endif lpencil_in(i_b2)=.true. lpencil_in(i_bb)=.true. if (delta_effect/=0.) lpencil_in(i_oxJ)=.true. if (meanfield_etat/=0.) then if (lweyl_gauge) then lpencil_in(i_jj)=.true. else lpencil_in(i_del2a)=.true. endif endif endif if (lpencil_in(i_del2A)) then if (lspherical_coords) then !WL: for the cylindrical case, lpencil_check says these pencils are not needed lpencil_in(i_jj)=.true. lpencil_in(i_graddivA)=.true. endif endif if (lpencil_in(i_uga)) then lpencil_in(i_aij)=.true. lpencil_in(i_uu)=.true. endif ! if (lpencil_in(i_ss12)) lpencil_in(i_sj)=.true. ! Pencils bij, del2a and graddiva come in a bundle. ! if ( lpencil_in(i_bij) .and. & ! (lpencil_in(i_del2a).or.lpencil_in(i_graddiva)) ) then ! lpencil_in(i_del2a)=.false. ! lpencil_in(i_graddiva)=.false. ! endif ! if (lpencil_in(i_del2a) .and. & ! (lpencil_in(i_bij).or.lpencil_in(i_graddiva)) ) then ! lpencil_in(i_bij)=.false. ! lpencil_in(i_graddiva)=.false. ! endif ! 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 ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! real, dimension (nx,3) :: bb_ext,bb_ext_pot real, dimension (nx) :: rho1_jxb,alpha_total real, dimension (nx) :: alpha_tmp real, dimension (nx) :: jcrossb2 real, dimension (nx) :: meanfield_Qs_func, meanfield_Qp_func real, dimension (nx) :: meanfield_Qs_der, meanfield_Qp_der, BiBk_Bki real, dimension (nx) :: meanfield_Bs21, meanfield_Bp21 real, dimension (nx) :: meanfield_Bs2, meanfield_Bp2 real, dimension (nx) :: meanfield_urms21,meanfield_etaB2 real, dimension (nx,3) :: Bk_Bki real :: B2_ext,c,s,kx 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,p%aa) ! bb if (lpencil(i_bb)) then call curl_mn(p%aij,p%bb,p%aa) ! ! 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.) then if (lbext_curvilinear.or.lcartesian_coords) 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.,1.,0.) is an azimuthal field in cylindrical ! coordinates and a polar one in spherical. ! if (omega_Bz_ext==0.) then B_ext_tmp=B_ext elseif (omega_Bz_ext/=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 else if (lcylindrical_coords) then if (omega_Bz_ext/=0.) & call fatal_error("calc_pencils_magnetic",& "precession of the external field not "//& "implemented for cylindrical coordinates") ! ! transform b_ext to other coordinate systems ! B_ext_tmp(1)= B_ext(1)*cos(y(m)) + B_ext(2)*sin(y(m)) B_ext_tmp(2)= -B_ext(1)*sin(y(m)) + B_ext(2)*cos(y(m)) B_ext_tmp(3)= B_ext(3) else if (lspherical_coords) then if (omega_Bz_ext/=0.) & call fatal_error("calc_pencils_magnetic",& "precession of the external field not "//& "implemented for spherical coordinates") B_ext_tmp(1)= B_ext(1)*sinth(m)*cos(z(n)) + B_ext(2)*sinth(m)*sin(z(n)) + B_ext(3)*costh(m) B_ext_tmp(2)= B_ext(1)*costh(m)*cos(z(n)) + B_ext(2)*costh(m)*sin(z(n)) - B_ext(3)*sinth(m) B_ext_tmp(3)=-B_ext(1) *sin(z(n)) + B_ext(2) *cos(z(n)) endif ! ! add the external field ! if (B_ext_tmp(1)/=0.) p%bb(:,1)=p%bb(:,1)+B_ext_tmp(1) if (B_ext_tmp(2)/=0.) p%bb(:,2)=p%bb(:,2)+B_ext_tmp(2) if (B_ext_tmp(3)/=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) ! uxb if (lpencil(i_uxb)) then call cross_mn(p%uu,p%bb,p%uxb) ! 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 ! uga ! DM : this requires later attention if (lpencil(i_uga)) then if (.not.lcartesian_coords) then call warning("calc_pencils_magnetic","u_dot_grad A not implemented for non-cartesian coordinates") else call u_dot_grad(f,iaa,p%aij,p%uu,p%uga,UPWIND=lupw_aa) endif endif ! ! 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 if (lcartesian_coords) then call gij_etc(f,iaa,p%aa,p%aij,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,p%bb) else call gij_etc(f,iaa,p%aa,p%aij,p%bij,GRADDIV=p%graddiva) call curl_mn(p%bij,p%jj,p%bb) if (lpencil(i_del2a)) p%del2a=p%graddiva-p%jj ! if (lpencil(i_del2a)) call del2v(f,iaa,p%del2a,p%aij,p%aa) endif 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 ! 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))<=tini).or.(abs(p%b2(ix))<=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) 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))<=tini).or.(abs(p%b2(ix))<=tini))then p%jperp=0 else p%jperp=sqrt(jcrossb2(ix))/sqrt(p%b2(ix)) endif enddo ! if (lpencil_check) then ! sinjb=sqrt(1-(modulo(p%cosjb + 1.0, 2.0) - 1)**2) ! else ! sinjb=sqrt(1-p%cosjb**2) ! endif ! p%jperp=sqrt(p%j2)*sinjb 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 if (lmeanfield_jxb) then if (lmeanfield_jxb_with_vA2) then meanfield_urms21=1./(3.*meanfield_kf*meanfield_etat)**2 meanfield_Qs_func=1.+(meanfield_Qs-1.)*(1.-2*pi_1*atan(p%vA2*meanfield_urms21)) meanfield_Qp_func=1.+(meanfield_Qp-1.)*(1.-2*pi_1*atan(p%vA2*meanfield_urms21)) meanfield_Qs_der=-2*pi_1*(meanfield_Qs-1.)/(1.+(p%vA2*meanfield_urms21)**2) meanfield_Qp_der=-2*pi_1*(meanfield_Qp-1.)/(1.+(p%vA2*meanfield_urms21)**2) call multsv_mn(meanfield_Qs_func,p%jxb,p%jxb) ! call multmv_transp(p%bij,p%bb,Bk_Bki) !call multsv_mn_add(meanfield_Qs_func-meanfield_Qp_func-p%b2*meanfield_Qp_der,Bk_Bki,p%jxb) !call multsv_mn_add(meanfield_Qs_func-meanfield_Qp_func-p%vA2*meanfield_urms21/meanfield_Bp2*meanfield_Qp_der,Bk_Bki,p%jxb) ! call multsv_mn_add(meanfield_Qs_func-meanfield_Qp_func,Bk_Bki,p%jxb) !call dot(Bk_Bki,p%bb,BiBk_Bki) ! call multsv_mn_add(2*meanfield_Qp_der*BiBk_Bki*p%rho1*meanfield_urms21,p%bb,p%jxb) else meanfield_Bs21=1./meanfield_Bs**2 meanfield_Bp21=1./meanfield_Bp**2 meanfield_Qs_func=1.+(meanfield_Qs-1.)*(1.-2*pi_1*atan(p%b2*meanfield_Bs21)) meanfield_Qp_func=1.+(meanfield_Qp-1.)*(1.-2*pi_1*atan(p%b2*meanfield_Bp21)) meanfield_Qs_der=-2*pi_1*(meanfield_Qs-1.)*meanfield_Bs21/(1.+(p%b2*meanfield_Bs21)**2) meanfield_Qp_der=-2*pi_1*(meanfield_Qp-1.)*meanfield_Bp21/(1.+(p%b2*meanfield_Bp21)**2) call multsv_mn(meanfield_Qs_func,p%jxb,p%jxb) call multmv_transp(p%bij,p%bb,Bk_Bki) call multsv_mn_add(meanfield_Qs_func-meanfield_Qp_func-p%b2*meanfield_Qp_der,Bk_Bki,p%jxb) call dot(Bk_Bki,p%bb,BiBk_Bki) call multsv_mn_add(2*meanfield_Qs_der*BiBk_Bki,p%bb,p%jxb) endif 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))<=tini).or.(abs(p%b2(ix))<=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) ! beta1 if (lpencil(i_beta1)) p%beta1=0.5*p%b2*mu01/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) ! del4a if (lpencil(i_del4a)) call del4v(f,iaa,p%del4a) ! del6a if (lpencil(i_del6a)) call del6v(f,iaa,p%del6a) ! 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)) ! ! mf_EMF ! needed if a mean field (mf) model is calculated ! if (lpencil(i_mf_EMF)) then kx=2*pi/Lx ! ! possibility of dynamical alpha ! if (lalpm.and..not.lmeanfield_noalpm) then alpha_total=alpha_effect*alpha_tmp+f(l1:l2,m,n,ialpm) else alpha_total=alpha_effect*alpha_tmp endif ! ! possibility of conventional alpha quenching (rescales alpha_total) ! initialize EMF with alpha_total*bb ! if (alpha_quenching/=0.) alpha_total=alpha_total/(1.+alpha_quenching*p%b2) call multsv_mn(alpha_total,p%bb,p%mf_EMF) ! ! add possible delta x J effect and turbulent diffusion to EMF ! if (delta_effect/=0.) p%mf_EMF=p%mf_EMF+delta_effect*p%oxJ if (meanfield_etat/=0.) then if (lweyl_gauge) then if (meanfield_etaB/=0.) then meanfield_etaB2=meanfield_etaB**2 call multsv_mn_add(meanfield_etat/sqrt(1.+p%b2/meanfield_etaB2),p%jj,p%mf_EMF) else p%mf_EMF=p%mf_EMF-meanfield_etat*p%jj endif else p%mf_EMF=p%mf_EMF+meanfield_etat*p%del2a endif ! ! allow for possibility of variable etat ! if (ietat/=0) then call multsv_mn_add(-f(l1:l2,m,n,ietat),p%jj,p%mf_EMF) endif endif endif if (lpencil(i_mf_EMFdotB)) call dot_mn(p%mf_EMF,p%bb,p%mf_EMFdotB) ! ! 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(f,df,p) ! ! magnetic field evolution ! ! calculate dA/dt=uxB+3/2 Omega_0 A_y x_dir -eta mu_0 J ! for mean field calculations one can also add dA/dt=...+alpha*bb+delta*WXJ ! add jxb/rho to momentum equation ! add eta mu_0 j2/rho to entropy equation ! ! 7-oct-09/axel: adapted from magnetic ! use Debug_IO, only: output_pencil use Deriv, only: der6 use Diagnostics use EquationOfState, only: eoscalc,gamma_m1 use Mpicomm, only: stop_it use Special, only: special_calc_magnetic use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx,3) :: bpol real, dimension (nx) :: aphi,bphi,d2aphi,d2bphi,bpol2 real, dimension (nx) :: alpha_tmp,alpha_tmp2,Omega_tmp real :: eta_tot integer :: i integer, parameter :: nxy=nxgrid*nygrid ! intent(inout) :: f,p intent(inout) :: df ! ! identify module and boundary conditions ! if (headtt.or.ldebug) print*,'daa_dt: SOLVE' if (headtt) then call identify_bcs('Aphi',iaphi) call identify_bcs('Bphi',ibphi) endif ! ! define aphi and bphi, and calculate bpol ! aphi=f(l1:l2,m,n,iaphi) bphi=f(l1:l2,m,n,ibphi) call curl_horizontal(f,iaphi,Bpol) ! ! Diffusion operator ! call del2(f,iaphi,d2aphi) call del2(f,ibphi,d2bphi) ! ! calculate D2a = del2a - 1/pomega^2 and same for b ! d2aphi=d2aphi-aphi*r2_mn*sin2th(m) d2bphi=d2bphi-bphi*r2_mn*sin2th(m) ! ! add diffusion term to the right-hand side. ! Use total eta, eta_tot, which includes microscopic and meanfield eta. ! eta_tot=eta+meanfield_etat df(l1:l2,m,n,iaphi)=df(l1:l2,m,n,iaphi)+eta_tot*d2aphi df(l1:l2,m,n,ibphi)=df(l1:l2,m,n,ibphi)+eta_tot*d2bphi ! ! add alpha effect, note that j=-D2a ! alpha_tmp=alpha_effect*alpha_x*alpha_y(m) alpha_tmp2=alpha_effect*(dalpha_x*Bpol(:,2)-r1_mn*dalpha_y(m)*Bpol(:,1)) ! ! Add alpha effect. At the moment this ignores the grad(alpha) terms ! df(l1:l2,m,n,iaphi)=df(l1:l2,m,n,iaphi)+alpha_tmp*bphi if (lalpha_phi_equation) then df(l1:l2,m,n,ibphi)=df(l1:l2,m,n,ibphi)-alpha_tmp*d2aphi+alpha_tmp2 endif ! ! differential rotation, need poloidal field for this, and ! add pomega*Bpol.grad(Omega) to the dBphi/dt equation. ! Note that grad_theta(Omega)=r1_mn*dOmega_y, but r1_mn cancels with r_mn. ! if (lOmega_effect) then !! Omega_tmp=sinth(m)*(r_mn*Bpol(:,1)*dOmega_x+Bpol(:,2)*dOmega_y(m)) Omega_tmp=sinth(m)*(r_mn*Bpol(:,1)*dOmega_x*Omega_y(m)+2*Bpol(:,2)*dOmega_y(m)*Omega_x) df(l1:l2,m,n,ibphi)=df(l1:l2,m,n,ibphi)+Omega_ampl*Omega_tmp endif ! ! allow for special routines ! if (lspecial) call special_calc_magnetic(f,df,p) ! ! Multiply resistivity by Nyquist scale, for resistive time-step. ! We include possible contribution from meanfield_etat, which is however ! only invoked in mean field models. ! Allow for variable etat (mean field theory) ! if (lfirst.and.ldt) then diffus_eta=diffus_eta+eta_tot diffus_eta=diffus_eta*dxyz_2 endif ! if (headtt.or.ldebug) then print*, 'daa_dt: max(diffus_eta) =', maxval(diffus_eta) endif ! ! if (linduction) & ! df(l1:l2,m,n,iax:iaz) = df(l1:l2,m,n,iax:iaz) + uxb_upw + fres ! endif ! ! Alpha effect ! additional terms if Mean Field Theory is included ! ! if (lmeanfield_theory.and. & ! (meanfield_etat/=0..or.alpha_effect/=0..or.delta_effect/=0.)) then ! df(l1:l2,m,n,iax:iaz)=df(l1:l2,m,n,iax:iaz)+p%mf_EMF ! if (lOmega_effect) call Omega_effect(f,df) ! endif ! if (ldiagnos) then if (idiag_aphi2m/=0) call sum_mn_name(aphi**2,idiag_aphi2m) if (idiag_bphi2m/=0) call sum_mn_name(bphi**2,idiag_bphi2m) if (idiag_bpol2m/=0) then call dot2(Bpol,Bpol2) call sum_mn_name(Bpol2,idiag_bpol2m) endif endif ! endsubroutine daa_dt !*********************************************************************** subroutine time_integrals_magnetic(f,p) ! ! Calculate time_integrals within each pencil (as long as each ! pencil case p still contains the current data). This routine ! is now being called at the end of equ. ! ! 28-jun-07/axel+mreinhard: coded ! 24-jun-08/axel: moved call to this routine to the individual pde routines ! 1-jul-08/axel: moved this part to magnetic ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(inout) :: f intent(in) :: p ! if (ibbt/=0) f(l1:l2,m,n,ibxt:ibzt)=f(l1:l2,m,n,ibxt:ibzt)+dt*p%bb if (ijjt/=0) f(l1:l2,m,n,ijxt:ijzt)=f(l1:l2,m,n,ijxt:ijzt)+dt*p%jj ! endsubroutine time_integrals_magnetic !*********************************************************************** subroutine df_diagnos_magnetic(df,p) ! ! calculate diagnostics that involves df ! Here we calculate and . ! The latter is calculated as - ! This is used in dynamo theory for checking the minimal tau approximation. ! ! 10-oct-06/axel: coded ! use Diagnostics use Sub ! real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx,3) :: uudot,aadot,udotxb,B1_gradu real, dimension (nx) :: B1dot_udotxb,B1dot_uxbdot,B1dot_aadot,uxbdot2 ! intent(in) :: df, p ! endsubroutine df_diagnos_magnetic !*********************************************************************** subroutine set_border_magnetic(f,df,p) ! ! Calculates the driving term for the border profile ! of the aa variable. ! ! 28-jul-06/wlad: coded ! use BorderProfiles, only: border_driving,set_border_initcond use Mpicomm, only: stop_it ! real, dimension(mx,my,mz,mfarray) :: f type (pencil_case) :: p real, dimension(mx,my,mz,mvar) :: df real, dimension(nx,3) :: f_target integer :: ju,j ! select case (borderaa) ! case ('zero','0') f_target=0. ! case ('initial-condition') do j=1,3 ju=j+iaa-1 call set_border_initcond(f,ju,f_target(:,j)) enddo ! case ('nothing') if (lroot.and.ip<=5) & print*,"set_border_magnetic: borderaa='nothing'" ! case default write(unit=errormsg,fmt=*) & 'set_border_magnetic: No such value for borderaa: ', & trim(borderaa) call fatal_error('set_border_magnetic',errormsg) endselect ! if (borderaa /= 'nothing') then do j=1,3 ju=j+iaa-1 call border_driving(f,df,p,f_target(:,j),ju) enddo endif ! endsubroutine set_border_magnetic !*********************************************************************** subroutine eta_shell(p,eta_mn,geta) ! ! 24-nov-03/dave: coded ! 23-jun-09/axel: generalized to lcylinder_in_a_box ! use Sub, only: step, der_step use Mpicomm, only: stop_it ! type (pencil_case) :: p real, dimension (nx) :: eta_mn real, dimension (nx) :: prof,eta_r real, dimension (nx,3) :: geta real :: d_int,d_ext ! eta_r=0. ! if (eta_int > 0.) then d_int = eta_int - eta else d_int = 0. endif if (eta_ext > 0.) then d_ext = eta_ext - eta else d_ext = 0. endif ! ! calculate steps in resistivity ! make this dependent on the geometry used ! ! (i) lcylinder_in_a_box ! if (lcylinder_in_a_box.or.lcylindrical_coords) then prof=step(p%rcyl_mn,r_int,wresistivity) eta_mn=d_int*(1-prof) prof=step(p%rcyl_mn,r_ext,wresistivity) eta_mn=eta+eta_mn+d_ext*prof ! ! calculate radial derivative of steps and gradient of eta ! prof=der_step(p%rcyl_mn,r_int,wresistivity) eta_r=-d_int*prof prof=der_step(p%rcyl_mn,r_ext,wresistivity) eta_r=eta_r+d_ext*prof geta=p%evr*spread(eta_r,2,3) ! ! (ii) lsphere_in_a_box ! elseif (lsphere_in_a_box.or.lspherical_coords) then prof=step(p%r_mn,r_int,wresistivity) eta_mn=d_int*(1-prof) prof=step(p%r_mn,r_ext,wresistivity) eta_mn=eta+eta_mn+d_ext*prof ! ! calculate radial derivative of steps and gradient of eta ! prof=der_step(p%r_mn,r_int,wresistivity) eta_r=-d_int*prof prof=der_step(p%r_mn,r_ext,wresistivity) eta_r=eta_r+d_ext*prof geta=p%evr*spread(eta_r,2,3) ! ! (iii) other cases are not implemented yet ! else call stop_it("eta_shell works only for spheres or cylinders") endif ! endsubroutine eta_shell !*********************************************************************** subroutine calc_bthresh() ! ! calculate bthresh from brms, give warnings if there are problems ! ! 6-aug-03/axel: coded ! ! give warning if brms is not set in prints.in ! if (idiag_brms==0) then if (lroot.and.lfirstpoint) then print*,'calc_bthresh: need to set brms in print.in to get bthresh' endif endif ! ! if nvec exceeds nbvecmax (=1/4) of points per processor, then begin to ! increase scaling factor on bthresh. These settings will stay in place ! until the next restart ! if (nbvec>nbvecmax.and.lfirstpoint) then print*,'calc_bthresh: processor ',iproc,': bthresh_scl,nbvec,nbvecmax=', & bthresh_scl,nbvec,nbvecmax bthresh_scl=bthresh_scl*1.2 endif ! ! calculate bthresh as a certain fraction of brms ! bthresh=bthresh_scl*bthresh_per_brms*brms ! endsubroutine calc_bthresh !*********************************************************************** subroutine rescaling_magnetic(f) ! ! Rescale magnetic field by factor rescale_aa, ! ! 22-feb-05/axel: coded ! 10-feb-09/petri: adapted from testfield ! use Sub, only: update_snaptime, read_snaptime ! real, dimension (mx,my,mz,mfarray) :: f character (len=fnlen) :: file logical :: lmagnetic_out logical, save :: lfirst_call=.true. ! intent(inout) :: f ! ! Reinitialize aa periodically if requested ! if (lreset_aa) then file=trim(datadir)//'/treset_aa.dat' if (lfirst_call) then call read_snaptime(trim(file),taareset,naareset,daareset,t) if (taareset==0 .or. taareset < t-daareset) then taareset=t+daareset endif lfirst_call=.false. endif ! ! Rescale when the time has come ! (Note that lmagnetic_out and ch are not used here) ! if (t >= taareset) then f(:,:,:,iax:iaz)=rescale_aa*f(:,:,:,iax:iaz) call update_snaptime(file,taareset,naareset,daareset,t,lmagnetic_out) endif endif ! endsubroutine rescaling_magnetic !*********************************************************************** subroutine calc_tau_aa_exterior(f,df) ! ! magnetic field relaxation to zero on time scale tau_aa_exterior within ! exterior region. For the time being this means z > zgrav. ! ! 29-jul-02/axel: coded ! use Gravity ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real :: scl integer :: j ! intent(in) :: f intent(out) :: df ! if (headtt) print*,'calc_tau_aa_exterior: tau=',tau_aa_exterior if (z(n)>zgrav) then scl=1./tau_aa_exterior do j=iax,iaz df(l1:l2,m,n,j)=df(l1:l2,m,n,j)-scl*f(l1:l2,m,n,j) enddo endif ! endsubroutine calc_tau_aa_exterior !*********************************************************************** subroutine read_magnetic_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=magnetic_init_pars, IOSTAT=iostat) ! 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(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=magnetic_run_pars, IOSTAT=iostat) ! 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(ix_loc,m1:m2 ,n1:n2 ,iax-1+slices%index) slices%xz =f(l1:l2 ,iy_loc,n1:n2 ,iax-1+slices%index) slices%xy =f(l1:l2 ,m1:m2 ,iz_loc ,iax-1+slices%index) slices%xy2=f(l1:l2 ,m1:m2 ,iz2_loc,iax-1+slices%index) if (lwrite_slice_xy3) & slices%xy3=f(l1:l2,m1:m2,iz3_loc,iax-1+slices%index) if (lwrite_slice_xy4) & slices%xy4=f(l1:l2,m1:m2,iz4_loc,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 (lwrite_slice_xy3) slices%xy3=>bb_xy3(:,:,slices%index) if (lwrite_slice_xy4) slices%xy4=>bb_xy4(:,:,slices%index) if (slices%index<=3) slices%ready=.true. endif ! ! Magnetic field (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 (lwrite_slice_xy3) slices%xy3=>jj_xy3(:,:,slices%index) if (lwrite_slice_xy4) slices%xy4=>jj_xy4(:,:,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 if (lwrite_slice_xy3) slices%xy3=>b2_xy3 if (lwrite_slice_xy4) slices%xy4=>b2_xy4 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 if (lwrite_slice_xy3) slices%xy3=>jb_xy3 if (lwrite_slice_xy4) slices%xy4=>jb_xy4 slices%ready=.true. ! ! Plasma beta ! case ('beta1') slices%yz =>beta1_yz slices%xz =>beta1_xz slices%xy =>beta1_xy slices%xy2=>beta1_xy2 if (lwrite_slice_xy3) slices%xy3=>beta1_xy3 if (lwrite_slice_xy4) slices%xy4=>beta1_xy4 slices%ready=.true. ! endselect ! endsubroutine get_slices_magnetic !*********************************************************************** subroutine calc_mfield ! ! calculate mean magnetic field from xy- or z-averages ! ! 19-jun-02/axel: moved from print to here ! 9-nov-02/axel: corrected bxmy(m,j); it used bzmy instead! ! use Mpicomm use Sub ! ! For vector output (of bb vectors) we need brms ! on all processors. It suffices to have this for times when lout=.true., ! but we need to broadcast the result to all procs. ! ! calculate brms (this requires that brms is set in print.in) ! broadcast result to other processors ! ! The following calculation involving spatial averages ! endsubroutine calc_mfield !*********************************************************************** subroutine alfven_x(ampl,f,iuu,iaa,ilnrho,kx) ! ! Alfven wave propagating in the x-direction ! ! ux = +sink(x-vA*t) ! Az = -cosk(x-vA*t)*sqrt(rho*mu0)/k ! ! Alfven and slow magnetosonic are the same here and both incompressible, and ! a fast magnetosonic (compressible) wave is also excited, but decoupled. ! ! satisfies the four equations ! dlnrho/dt = -ux' ! dux/dt = -cs2*(lnrho)' ! duy/dt = B0*By' ==> duy/dt = -B0*Az'' ! dBy/dt = B0*uy' ==> dAz/dt = -B0*ux ! ! 8-nov-03/axel: coded ! 29-apr-03/axel: added sqrt(rho*mu0)/k factor ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (nx) :: rho,ampl_Az real :: ampl,kx,ampl_lr,ampl_ux,ampl_uy integer :: iuu,iaa,ilnrho ! ! Amplitude factors ! ampl_lr=+0. ampl_ux=+0. ampl_uy=+ampl ! ! ux and Ay. ! Don't overwrite the density, just add to the log of it. ! do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,ilnrho)=ampl_lr*(sin(kx*x(l1:l2))+f(l1:l2,m,n,ilnrho)) f(l1:l2,m,n,iuu+0 )=ampl_ux*sin(kx*x(l1:l2)) f(l1:l2,m,n,iuu+1 )=ampl_uy*sin(kx*x(l1:l2)) rho=exp(f(l1:l2,m,n,ilnrho)) ampl_Az=-ampl*sqrt(rho*mu0)/kx f(l1:l2,m,n,iaa+2 )=ampl_Az*cos(kx*x(l1:l2)) enddo; enddo ! endsubroutine alfven_x !*********************************************************************** subroutine alfven_y(ampl,f,iuu,iaa,ky,mu0) ! ! Alfven wave propagating in the y-direction; can be used in 2-d runs. ! ux = cos(ky-ot), for B0y=1 and rho=1. ! Az = sin(ky-ot), ie Bx=-cos(ky-ot) ! ! [wd nov-2006: There should be a 1/ky in the aa term here and in ! alfven_x, I think] ! ! satisfies the equations ! dux/dt = Bx' ==> dux/dt = -Az'' ! dBx/dt = ux' ==> dAz/dt = -ux. ! ! 06-dec-06/wolf: adapted from alfven_z ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,ky,mu0 integer :: iuu,iaa ! ! ux and Az ! do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+0) = +ampl*cos(ky*y(m)) f(l1:l2,m,n,iaa+2) = -ampl*sin(ky*y(m))*sqrt(mu0)/ky enddo; enddo ! endsubroutine alfven_y !*********************************************************************** subroutine alfven_z(ampl,f,iuu,iaa,kz,mu0) ! ! Alfven wave propagating in the z-direction ! ux = cos(kz-ot), for B0z=1 and rho=1. ! Ay = sin(kz-ot), ie Bx=-cos(kz-ot) ! ! satisfies the equations ! dux/dt = Bx' ==> dux/dt = -Ay'' ! dBx/dt = ux' ==> dAy/dt = -ux. ! ! 18-aug-02/axel: coded ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,kz,mu0 integer :: iuu,iaa ! ! ux and Ay ! do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+0)=+ampl*cos(kz*z(n)) f(l1:l2,m,n,iaa+1)=+ampl*sin(kz*z(n))*sqrt(mu0) enddo; enddo ! endsubroutine alfven_z !*********************************************************************** subroutine alfven_xy(ampl,f,iuu,iaa,kx,ky) ! ! Alfven wave propagating in the xy-direction; can be used in 2-d runs. ! uz = cos(kx*x+ky*y-ot), for B0=(1,1,0) and rho=1. ! Ax = sin(kx*x+ky*y-ot), ! Ay = sin(kx*x+ky*y-ot), ! ! 16-jun-07/axel: adapted from alfven_y ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,kx,ky,om real, parameter :: mu0=1. integer :: iuu,iaa ! ! set ux, Ax, and Ay ! om=B_ext(1)*kx+B_ext(2)*ky do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+2)=+ampl*cos(kx*x(l1:l2)+ky*y(m)) f(l1:l2,m,n,iaa+0)=+ampl*sin(kx*x(l1:l2)+ky*y(m))*sqrt(mu0)/om*B_ext(2) f(l1:l2,m,n,iaa+1)=-ampl*sin(kx*x(l1:l2)+ky*y(m))*sqrt(mu0)/om*B_ext(1) enddo; enddo ! endsubroutine alfven_xy !*********************************************************************** subroutine alfven_xz(ampl,f,iuu,iaa,kx,kz) ! ! Alfven wave propagating in the xz-direction; can be used in 2-d runs. ! uz = cos(kx*x+kz*z-ot), for B0=(1,1,0) and rho=1. ! Ax = sin(kx*x+kz*z-ot), ! Az = sin(kx*x+kz*z-ot), ! ! 16-jun-07/axel: adapted from alfven_xy ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,kx,kz,om real, parameter :: mu0=1. integer :: iuu,iaa ! ! set ux, Ax, and Az ! om=B_ext(1)*kx+B_ext(3)*kz do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+2)=+ampl*cos(kx*x(l1:l2)+kz*z(n)) f(l1:l2,m,n,iaa+0)=+ampl*sin(kx*x(l1:l2)+kz*z(n))*sqrt(mu0)/om*B_ext(2) f(l1:l2,m,n,iaa+2)=-ampl*sin(kx*x(l1:l2)+kz*z(n))*sqrt(mu0)/om*B_ext(1) enddo; enddo ! endsubroutine alfven_xz !*********************************************************************** subroutine alfvenz_rot(ampl,f,iuu,iaa,kz,O) ! ! Alfven wave propagating in the z-direction (with Coriolis force) ! ux = cos(kz-ot), for B0z=1 and rho=1. ! Ay = sin(kz-ot), ie Bx=-cos(kz-ot) ! ! satisfies the equations ! dux/dt - 2Omega*uy = -Ay'' ! duy/dt + 2Omega*ux = +Ax'' ! dAx/dt = +uy ! dAy/dt = -ux ! ! 18-aug-02/axel: coded ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,kz,O,fac integer :: iuu,iaa ! ! ux, uy, Ax and Ay ! if (lroot) print*,'alfvenz_rot: Alfven wave with rotation; O,kz=',O,kz fac=-O+sqrt(O**2+kz**2) do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+0)=-ampl*sin(kz*z(n))*fac/kz f(l1:l2,m,n,iuu+1)=-ampl*cos(kz*z(n))*fac/kz f(l1:l2,m,n,iaa+0)=+ampl*sin(kz*z(n))/kz f(l1:l2,m,n,iaa+1)=+ampl*cos(kz*z(n))/kz enddo; enddo ! endsubroutine alfvenz_rot !*********************************************************************** subroutine alfvenz_rot_shear(ampl,f,iuu,iaa,kz,OO) ! ! Alfven wave propagating in the z-direction (with Coriolis force and shear) ! ! satisfies the equations ! dux/dt - 2*Omega*uy = -Ay'' ! duy/dt + (2-q)*Omega*ux = +Ax'' ! dAx/dt = q*Omega*Ay + uy ! dAy/dt = -ux ! ! Assume B0=rho0=mu0=1 ! ! 28-june-04/anders: coded ! real, dimension (mx,my,mz,mfarray) :: f real :: ampl,kz,OO complex :: fac integer :: iuu,iaa ! ! ux, uy, Ax and Ay ! if (lroot) print*,'alfvenz_rot_shear: '// & 'Alfven wave with rotation and shear; OO,kz=',OO,kz fac=cmplx(OO-sqrt(16*kz**2+OO**2),0.) do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,iuu+0)=f(l1:l2,m,n,iuu+0)+ampl*fac/(4*kz)*sin(kz*z(n)) f(l1:l2,m,n,iuu+1)=f(l1:l2,m,n,iuu+1)+ampl*real(exp(cmplx(0,z(n)*kz))* & fac*sqrt(2*kz**2+OO*fac)/(sqrt(2.)*kz*(-6*OO-fac))) f(l1:l2,m,n,iaa+0)=ampl*sin(kz*z(n))/kz f(l1:l2,m,n,iaa+1)=-ampl*2*sqrt(2.)*aimag(exp(cmplx(0,z(n)*kz))* & sqrt(2*kz**2+OO*fac)/(-6*OO-fac)/(cmplx(0,kz))) enddo; enddo ! endsubroutine alfvenz_rot_shear !*********************************************************************** subroutine fluxrings(ampl,f,ivar1,ivar2,profile) ! ! Magnetic flux rings. Constructed from a canonical ring which is the ! rotated and translated: ! AA(xxx) = D*AA0(D^(-1)*(xxx-xxx_disp)) , ! where AA0(xxx) is the canonical ring and D the rotation matrix ! corresponding to a rotation by phi around z, followed by a ! rotation by theta around y. ! The array was already initialized to zero before calling this ! routine. ! Optional argument `profile' allows to choose a different profile (see ! norm_ring()) ! use Mpicomm, only: stop_it ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (nx,3) :: tmpv real, dimension (nx) :: xx1,yy1,zz1 real, dimension(3) :: axis,disp real :: ampl,phi,theta,ct,st,cp,sp real :: fring,Iring,R0,width integer :: i,ivar,ivar1,ivar2,ivar3 character (len=*), optional :: profile character (len=labellen) :: prof ! if (present(profile)) then prof = profile else prof = 'tanh' endif ! ! fix ivar3=ivar1 (for now) ! ivar3=ivar1 ! ! initialize each ring ! if (any((/fring1,fring2,Iring1,Iring2/) /= 0.)) then ! fringX is the magnetic flux, IringX the current if (lroot) then print*, 'fluxrings: Initialising magnetic flux rings' endif do i=1,nrings if (i==1) then fring = fring1 ! magnetic flux along ring Iring = Iring1 ! current along ring (for twisted flux tube) R0 = Rring1 ! radius of ring width = wr1 ! ring thickness axis = axisr1 ! orientation disp = dispr1 ! position ivar = ivar1 elseif (i==2) then fring = fring2 Iring = Iring2 R0 = Rring2 width = wr2 axis = axisr2 disp = dispr2 ivar = ivar2 elseif (i==3) then fring = fring3 Iring = Iring3 R0 = Rring3 width = wr3 axis = axisr3 disp = dispr3 ivar = ivar3 else call stop_it('fluxrings: nrings is too big') endif phi = atan2(axis(2),axis(1)+epsi) theta = atan2(sqrt(axis(1)**2+axis(2)**2)+epsi,axis(3)) ct = cos(theta); st = sin(theta) cp = cos(phi) ; sp = sin(phi) ! Calculate D^(-1)*(xxx-disp) do n=n1,n2; do m=m1,m2 xx1= ct*cp*(x(l1:l2)-disp(1))+ct*sp*(y(m)-disp(2))-st*(z(n)-disp(3)) yy1=- sp*(x(l1:l2)-disp(1))+ cp*(y(m)-disp(2)) zz1= st*cp*(x(l1:l2)-disp(1))+st*sp*(y(m)-disp(2))+ct*(z(n)-disp(3)) call norm_ring(xx1,yy1,zz1,fring,Iring,R0,width,tmpv,PROFILE=prof) ! calculate D*tmpv f(l1:l2,m,n,ivar ) = f(l1:l2,m,n,ivar ) + ampl*( & + ct*cp*tmpv(:,1) - sp*tmpv(:,2) + st*cp*tmpv(:,3)) f(l1:l2,m,n,ivar+1) = f(l1:l2,m,n,ivar+1) + ampl*( & + ct*sp*tmpv(:,1) + cp*tmpv(:,2) + st*sp*tmpv(:,3)) f(l1:l2,m,n,ivar+2) = f(l1:l2,m,n,ivar+2) + ampl*( & - st *tmpv(:,1) + ct *tmpv(:,3)) enddo; enddo enddo endif if (lroot) print*, 'fluxrings: Magnetic flux rings initialized' ! endsubroutine fluxrings !*********************************************************************** subroutine norm_ring(xx1,yy1,zz1,fring,Iring,R0,width,vv,profile) ! ! Generate vector potential for a flux ring of magnetic flux FRING, ! current Iring (not correctly normalized), radius R0 and thickness ! WIDTH in normal orientation (lying in the x-y plane, centred at (0,0,0)). ! ! 1-may-02/wolf: coded (see the manual, Section C.3) ! 7-jun-09/axel: added gaussian and constant (or box) profiles ! use Mpicomm, only: stop_it use Sub ! real, dimension (nx,3) :: vv real, dimension (nx) :: xx1,yy1,zz1,phi,tmp real :: fring,Iring,R0,width character (len=*) :: profile ! ! magnetic ring, define r-R ! tmp = sqrt(xx1**2+yy1**2)-R0 ! ! choice of different profile functions ! select case (profile) ! ! gaussian profile, exp(-.5*(x/w)^2)/(sqrt(2*pi)*eps), ! so its derivative is .5*(1.+erf(-x/(sqrt(2)*eps)) ! case ('gaussian') vv(:,3) = - fring * .5*(1.+erfunc(tmp/(sqrt(2.)*width))) & * exp(-.5*(zz1/width)**2)/(sqrt(2.*pi)*width) ! ! tanh profile, so the delta function is approximated by 1/cosh^2. ! The name tanh is misleading, because the actual B frofile is ! 1./cosh^2, but this is harder to write. ! case ('tanh') vv(:,3) = - fring * 0.5*(1+tanh(tmp/width)) & * 0.5/width/cosh(zz1/width)**2 ! ! constant profile, so the delta function is approximated by the function ! delta(x) = 1/2w, if -w < x < w. ! case ('const') vv(:,3) = - fring * 0.5*(1.+max(-1.,min(tmp/width,1.))) & * 0.25/width*(1.-sign(1.,abs(zz1)-width)) ! ! there is no default option here ! case default call stop_it('norm_ring: No such fluxtube profile') endselect ! ! current ring (to twist the B-lines) ! tmp = width - sqrt(tmp**2 + zz1**2) tmp = Iring*0.5*(1+tanh(tmp/width)) ! Now the A_phi component phi = atan2(yy1,xx1) vv(:,1) = - tmp*sin(phi) vv(:,2) = tmp*cos(phi) ! endsubroutine norm_ring !*********************************************************************** subroutine torus_test(ampl,f) ! ! Initial field concentrated along torus well inside the computational ! domain. ! Implements the same field for cartesian and spherical cordinates. ! The field is of mixed parity (bb_pol symmetric, bb_tor antisymmetric) ! and the relative contributions of toroidal and poloidal field are ! determined by ! ampl(1) -- bb_pol (through aa_tor) ! ampl(3) -- bb_tor (through aa_pol) ! Uses x_max as reference radius. ! ! 05-may-2008/wolf: coded ! real :: ampl real, dimension (mx,my,mz,mfarray) :: f ! real, dimension (nx) :: xxi2,ee real, dimension (nx) :: costh,sinth,cosphi,sinphi,ss,rr,aar,aap real :: radius,width,r_cent ! intent(in) :: ampl intent(inout) :: f ! radius = xyz1(1) width = 0.1*radius r_cent = 0.6*radius ! if (lspherical_coords) then do n=n1,n2; do m=m1,m2 xxi2 = (x(l1:l2)*sin(y(m))-r_cent)**2 + x(l1:l2)**2*cos(y(m))**2 ee = ampl * exp(-0.5 * xxi2 / width**2) f(l1:l2,m,n,iax) = f(l1:l2,m,n,iax) + ee * x(l1:l2)*cos(y(m)) f(l1:l2,m,n,iaz) = f(l1:l2,m,n,iaz) + ee enddo; enddo else do n=n1,n2; do m=m1,m2 xxi2 = (sqrt(x(l1:l2)**2+y(m)**2) - r_cent)**2 + z(n)**2 ee = ampl * exp(-0.5 * xxi2 / width**2) aar = z(n) * ee aap = ee ss = sqrt(x(l1:l2)**2+y(m)**2) rr = sqrt(x(l1:l2)**2+y(m)**2+z(n)**2) ss = max(ss, tini) rr = max(rr, tini) costh = z(n)/rr sinth = ss/rr cosphi = x(l1:l2)/ss sinphi = y(m)/ss f(l1:l2,m,n,iax) = f(l1:l2,m,n,iax) + aar*sinth*cosphi - aap*sinphi f(l1:l2,m,n,iay) = f(l1:l2,m,n,iay) + aar*sinth*sinphi + aap*cosphi f(l1:l2,m,n,iaz) = f(l1:l2,m,n,iaz) + aar*costh enddo; enddo endif ! endsubroutine torus_test !*********************************************************************** subroutine force_free_jet(mu) ! ! Force free magnetic field configuration for jet simulations ! with a fixed accretion disc at the bottom boundary. ! ! The input parameter mu specifies the radial dependency of ! the magnetic field in the disc. ! ! Solves the laplace equation in cylindrical coordinates for the ! phi-component of the vector potential. A_r and A_z are taken to ! be zero. ! ! nabla**2 A_phi - A_phi / r**2 = 0 ! ! For the desired boundary condition in the accretion disc ! ! B_r=B0*r**(mu-1) (z == 0) ! ! the solution is ! ! A_phi = Hypergeometric2F1( (1-mu)/2, (2+mu)/2, 2, xi**2 ) ! *xi*(r**2+z**2)**(mu/2) ! ! where xi = sqrt(r**2/(r**2+z**2)) ! ! ! 30-may-04/tobi: coded ! use Sub, only: hypergeometric2F1,gamma_function use Deriv, only: der use Debug_IO, only: output ! real, intent(in) :: mu real :: xi2,A_phi real :: r2 real :: B1r_,B1z_,B1 real, parameter :: tol=10*epsilon(1.0) integer :: l real, dimension(mx,my,mz) :: Ax_ext,Ay_ext real, dimension(nx,3) :: bb_ext_pot real, dimension(nx) :: bb_x,bb_y,bb_z ! ! calculate un-normalized |B| at r=r_ref and z=0 for later normalization ! if (lroot.and.ip<=5) print*,'FORCE_FREE_JET: calculating normalization' ! B1r_=sin(pi*mu/2)*gamma_function( abs(mu) /2) / & gamma_function((1+abs(mu))/2) ! B1z_=cos(pi*mu/2)*gamma_function((1+abs(mu))/2) / & gamma_function((2+abs(mu))/2) ! B1=sqrt(4/pi)*r_ref**(mu-1)*sqrt(B1r_**2+B1z_**2) ! ! calculate external vector potential ! if (lroot) print*,'FORCE_FREE_JET: calculating external vector potential' ! if (lforce_free_test) then ! if (lroot) print*,'FORCE_FREE_JET: using analytic solution for mu=-1' do l=1,mx; do m=1,my; do n=1,mz Ax_ext=-2*y(m)*(1-z(n)/sqrt(x(l)**2+y(m)**2+z(n)**2))/(x(l)**2+y(m)**2)/B1 Ay_ext= 2*x(l)*(1-z(n)/sqrt(x(l)**2+y(m)**2+z(n)**2))/(x(l)**2+y(m)**2)/B1 enddo; enddo; enddo ! else ! do l=1,mx; do m=1,my; do n=1,mz r2=x(l)**2+y(m)**2 xi2=r2/(r2+z(n)**2) A_phi=hypergeometric2F1((1-mu)/2,(2+mu)/2,2.0,xi2,tol) & *sqrt(xi2)*sqrt(r2+z(n)**2)**mu/B1 ! Ax_ext(l,m,n)=-y(m)*A_phi/sqrt(r2) Ay_ext(l,m,n)= x(l)*A_phi/sqrt(r2) enddo; enddo; enddo ! endif ! ! calculate external magnetic field ! if (lroot.and.ip<=5) & print*,'FORCE_FREE_JET: calculating the external magnetic field' ! do n=n1,n2 do m=m1,m2 ! call der(Ay_ext,bb_x,3) ! bb_ext_pot(:,1)=-bb_x ! call der(Ax_ext,bb_y,3) ! bb_ext_pot(:,2)= bb_y ! call der(Ay_ext,bb_z,1) ! bb_ext_pot(:,3)= bb_z ! call der(Ax_ext,bb_z,2) ! bb_ext_pot(:,3)=bb_ext_pot(:,3)-bb_z ! call set_global(bb_ext_pot,m,n,'B_ext_pot',nx) enddo enddo ! if (ip<=5) then call output(trim(directory)//'/Ax_ext.dat',Ax_ext,1) call output(trim(directory)//'/Ay_ext.dat',Ay_ext,1) endif ! endsubroutine force_free_jet !*********************************************************************** subroutine piecew_dipole_aa(ampl,inclaa,f,ivar) ! ! A field that is vertical uniform for rR_ext, and potential in the shell R_intr_ext) ! r > r_ext sigma0 = a(0,3)*r_mn + b(0,3)*r_2_mn sigma1 = a(1,3)*r_mn + b(1,3)*r_2_mn endwhere where(r_mn= r_int) ar=ampl_B0*80.*2.*(3.*sin(theta_mn)**2-2.)* & ( 1./36.*r_mn**5 - 1./12.*(r_int+r_ext)*r_mn**4 & + 1./14.*(r_int**2+4.*r_int*r_ext+r_ext**2)*r_mn**3 & - 1./3.*(r_int**2*r_ext+r_int*r_ext**2)*r_mn**2 & - 1./25.*r_int**2*r_ext**2*r_mn & + 1./5.*r_int**2*r_ext**2*r_mn*log(r_mn) ) atheta=-ampl_B0*80.*sin(2.*theta_mn)* & ( 7./36.*r_mn**5 - 1./2.*(r_int+r_ext)*r_mn**4 & + 5./14.*(r_int**2+4.*r_int*r_ext+r_ext**2)*r_mn**3 & - 4./3.*(r_int**2*r_ext+r_int*r_ext**2)*r_mn**2 & + 2./25.*r_int**2*r_ext**2*r_mn & + 3./5.*r_int**2*r_ext**2*r_mn*log(r_mn) ) aphi=ampl_B0*5./8.*sin(theta_mn)* & ( 4.*r_ext*r_mn - 3.*r_mn**2 - r_int**4/r_mn**2 ) endwhere ! where (r_mn > r_ext) ar=C_ext*ampl_B0*80.*2.*(3.*sin(theta_mn)**2-2.)/r_mn**4 atheta=-2.*C_ext*ampl_B0*80.*sin(2.*theta_mn)/r_mn**4 aphi=ampl_B0*A_ext/r_mn**2*sin(theta_mn) endwhere ! ! debug checks -- look at a pencil near the centre... if (ip<=4 .and. imn==(ny+1)*nz/2) then print*,'r_int,r_ext',r_int,r_ext write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(r_mn), imn, iproc:', & iproc, imn, minval(r_mn), maxval(r_mn) write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(theta_mn), imn, iproc:', & iproc, imn, minval(theta_mn), maxval(theta_mn) write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(phi_mn), imn, iproc:', & iproc, imn, minval(phi_mn), maxval(phi_mn) write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(ar), imn, iproc:', & iproc, imn, minval(ar), maxval(ar) write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(atheta), imn, iproc:', & iproc, imn, minval(atheta), maxval(atheta) write(*,'(a45,2i6,2f15.7)') & 'geo_benchmark_B: minmax(aphi), imn, iproc:', & iproc, imn, minval(aphi), maxval(aphi) endif ! case ('geo-benchmark-case2') if (lroot .and. imn==1) print*, 'geo_benchmark_B: geo-benchmark-case2 not yet coded.' ! case default if (lroot .and. imn==1) print*,'geo_benchmark_B: case not defined!' call stop_it("") endselect enddo f(l1:l2,m,n,iax)=sin(theta_mn)*cos(phi_mn)*ar + cos(theta_mn)*cos(phi_mn)*atheta - sin(phi_mn)*aphi f(l1:l2,m,n,iay)=sin(theta_mn)*sin(phi_mn)*ar + cos(theta_mn)*sin(phi_mn)*atheta + cos(phi_mn)*aphi f(l1:l2,m,n,iaz)=cos(theta_mn)*ar - sin(theta_mn)*atheta enddo ! if (ip<=14) then print*,'geo_benchmark_B: minmax(ax) on iproc:', iproc, minval(f(l1:l2,m1:m2,n1:n2,iax)),maxval(f(l1:l2,m1:m2,n1:n2,iax)) print*,'geo_benchmark_B: minmax(ay) on iproc:', iproc, minval(f(l1:l2,m1:m2,n1:n2,iay)),maxval(f(l1:l2,m1:m2,n1:n2,iay)) print*,'geo_benchmark_B: minmax(az) on iproc:', iproc, minval(f(l1:l2,m1:m2,n1:n2,iaz)),maxval(f(l1:l2,m1:m2,n1:n2,iaz)) endif ! endsubroutine geo_benchmark_B !*********************************************************************** subroutine eta_rdep(eta_r,geta_r,rdep_profile) ! ! 2-jul-2009/koen: creates an r-dependent resistivity for RFP studies ! real, dimension(nx,my) :: eta_r,r2,rmax2,gradr_eta_r real, dimension(nx,my,3) :: geta_r character (len=labellen) :: rdep_profile integer :: i,j ! intent(out) :: eta_r,geta_r ! ! adapted from etazdep ! Note the usage of mixed array lengths (nx and my) ! select case (rdep_profile) case ('schnack89') do i=1,nx do j=1,my r2(i,j)=x(i+l1-1)**2+y(j)**2 rmax2=1. !value should come from start.in? enddo enddo ! ! define eta_r: resistivity profile from Y.L. Ho, S.C. Prager & ! D.D. Schnack, Phys rev letters vol 62 nr 13 1989 ! and define gradr_eta_r: 1/r *d_r(eta_r)) ! eta_r = eta*(1+9*(r2/rmax2)**15)**2 gradr_eta_r= 540*eta*(1+9*(r2/rmax2)**15)*(r2/rmax2)**14/rmax2**0.5 ! ! gradient ! do i=1,nx geta_r(i,:,1) = x(i+l1-1)*gradr_eta_r(i,:) enddo do i=1,my geta_r(:,i,2) = y(i)*gradr_eta_r(:,i) enddo geta_r(:,:,3) = 0. ! ! possibility of constant eta_r (as a test) ! case ('constant') eta_r=eta ! ! gradient geta_r(:,:,1) = 0. geta_r(:,:,2) = 0. geta_r(:,:,3) = 0. endselect ! endsubroutine eta_rdep !*********************************************************************** subroutine eta_zdep(eta_z,geta_z,zdep_profile) ! ! creates a z-dependent resistivity for protoplanetary disk studies ! ! 12-jul-2005/joishi: coded ! use General, only:erfcc ! real, dimension(mz) :: eta_z,z2 real, dimension(mz,3) :: geta_z character (len=labellen) :: zdep_profile ! integer :: i ! intent(out) :: eta_z,geta_z ! select case (zdep_profile) case ('fs') z2 = z**2. ! resistivity profile from Fleming & Stone (ApJ 585:908-920) eta_z = eta*exp(-z2/2.+sigma_ratio*erfcc(abs(z))/4.) ! ! its gradient: geta_z(:,1) = 0. geta_z(:,2) = 0. geta_z(:,3) = eta_z*(-z-sign(1.,z)*sigma_ratio*exp(-z2)/(2.*sqrt(pi))) ! case ('tanh') ! default to spread gradient over ~5 grid cells. if (eta_width == 0.) eta_width = 5.*dz eta_z = eta*0.5*(tanh((z + eta_z0)/eta_width) & - tanh((z - eta_z0)/eta_width)) ! ! its gradient: geta_z(:,1) = 0. geta_z(:,2) = 0. geta_z(:,3) = -eta/(2.*eta_width) * ((tanh((z + eta_z0)/eta_width))**2. & - (tanh((z - eta_z0)/eta_width))**2.) ! endselect ! endsubroutine eta_zdep !************************************************************************ subroutine bb_unitvec_shock(f,bb_hat) ! ! Compute unit vector along the magnetic field. ! Accurate to 2nd order. ! Tries to avoid division by zero. ! Taken from http://nuclear.llnl.gov/CNP/apt/apt/aptvunb.html. ! If anybody knows a more accurate way of doing this, please modify. ! ! 16-aug-06/tobi: coded ! real, dimension (mx,my,mz,mfarray), intent (in) :: f real, dimension (mx,3), intent (out) :: bb_hat ! !Tobi: Not sure about this value real, parameter :: tol=1e-11 ! real, dimension (mx,3) :: bb,bb2 real, dimension (mx) :: bb_len,aerr2 real :: fac integer :: j ! ! Compute magnetic field from vector potential. ! bb=0. ! if (nxgrid/=1) then fac = 1/(2*dx) bb(l1-2:l2+2,3) = bb(l1-2:l2+2,3) + fac*( f(l1-1:l2+3,m ,n ,iay) & - f(l1-3:l2+1,m ,n ,iay) ) bb(l1-2:l2+2,2) = bb(l1-2:l2+2,2) - fac*( f(l1-1:l2+3,m ,n ,iaz) & - f(l1-3:l2+1,m ,n ,iaz) ) endif ! if (nygrid/=1) then fac = 1/(2*dy) bb(l1-2:l2+2,1) = bb(l1-2:l2+2,1) + fac*( f(l1-2:l2+2,m+1,n ,iaz) & - f(l1-2:l2+2,m-1,n ,iaz) ) bb(l1-2:l2+2,3) = bb(l1-2:l2+2,3) - fac*( f(l1-2:l2+2,m+1,n ,iax) & - f(l1-2:l2+2,m-1,n ,iax) ) endif ! if (nzgrid/=1) then fac = 1/(2*dz) bb(l1-2:l2+2,2) = bb(l1-2:l2+2,2) + fac*( f(l1-2:l2+2,m ,n+1,iax) & - f(l1-2:l2+2,m ,n-1,iax) ) bb(l1-2:l2+2,1) = bb(l1-2:l2+2,1) - fac*( f(l1-2:l2+2,m ,n+1,iay) & - f(l1-2:l2+2,m ,n-1,iay) ) endif ! ! Add external magnetic field. ! do j=1,3; bb(:,j) = bb(:,j) + B_ext(j); enddo ! ! Truncate small components to zero. ! bb2 = bb**2 ! aerr2 = tol**2 * max(sum(bb2,2),1.) ! do j=1,3 where (bb2(:,j) < aerr2) bb_hat(:,j) = 0. elsewhere bb_hat(:,j) = bb(:,j) endwhere enddo ! ! Get unit vector. ! bb_len = sqrt(sum(bb_hat**2,2)) ! do j=1,3; bb_hat(:,j) = bb_hat(:,j)/(bb_len+tini); enddo ! endsubroutine bb_unitvec_shock !*********************************************************************** subroutine input_persistent_magnetic(id,lun,done) ! ! Read in the stored phase and amplitude for the ! correction of the Beltrami wave forcing ! ! 5-apr-08/axel: adapted from input_persistent_forcing ! 01-Jun-2015/Bourdin.KIS: reworked ! integer :: id,lun logical :: done ! select case (id) case (id_record_MAGNETIC_PHASE) if (read_persist ('MAGNETIC_PHASE', phase_beltrami)) return done = .true. case (id_record_MAGNETIC_AMPL) if (read_persist ('MAGNETIC_AMPL', ampl_beltrami)) return done = .true. endselect ! endsubroutine input_persistent_magnetic !*********************************************************************** logical function output_persistent_magnetic(lun) ! ! Write the stored phase and amplitude for the ! correction of the Beltrami wave forcing ! ! 5-apr-08/axel: adapted from output_persistent_forcing ! 01-Jun-2015/Bourdin.KIS: reworked ! use IO, only: write_persist ! integer :: lun ! if (lroot .and. (ip < 14)) then if (phase_beltrami >= 0.) & print *,'output_persistent_magnetic: ',phase_beltrami,ampl_beltrami endif ! ! write details ! output_persistent_magnetic = .true. ! if (write_persist ('MAGNETIC_PHASE', id_record_MAGNETIC_PHASE, phase_beltrami)) return if (write_persist ('MAGNETIC_AMPL', id_record_MAGNETIC_AMPL, ampl_beltrami)) return ! output_persistent_magnetic = .false. ! endfunction output_persistent_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 use FArrayManager, only: farray_index_append ! integer :: iname,inamex,inamey,inamez,ixy,ixz,irz,inamer,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_aphi2m=0; idiag_bphi2m=0; idiag_bpol2m=0 endif ! ! check for those quantities that we want to evaluate online ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'aphi2m',idiag_aphi2m) call parse_name(iname,cname(iname),cform(iname),'bphi2m',idiag_bphi2m) call parse_name(iname,cname(iname),cform(iname),'bpol2m',idiag_bpol2m) enddo ! ! write column, idiag_XYZ, where our variable XYZ is stored ! if (lwr) then call farray_index_append('nname',nname) call farray_index_append('nnamexy',nnamexy) call farray_index_append('nnamexz',nnamexz) call farray_index_append('nnamez',nnamez) call farray_index_append('iaphi',iaphi) call farray_index_append('ibphi',ibphi) endif ! endsubroutine rprint_magnetic !*********************************************************************** subroutine split_update_magnetic(f) ! ! dummy ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine split_update_magnetic !*********************************************************************** subroutine expand_shands_magnetic() ! ! Presently dummy, for later use ! endsubroutine expand_shands_magnetic !*********************************************************************** endmodule Magnetic