MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_mhd_phys.t
Go to the documentation of this file.
1 !> Magneto-hydrodynamics module
3 
4 #include "amrvac.h"
5 
6  use mod_global_parameters, only: std_len, const_c
10  use mod_physics
11  use mod_comm_lib, only: mpistop
13 
14  implicit none
15  private
16 
17  !> Whether an energy equation is used
18  logical, public, protected :: mhd_energy = .true.
19 
20  !> Whether thermal conduction is used
21  logical, public, protected :: mhd_thermal_conduction = .false.
22  !> type of fluid for thermal conduction
23  type(tc_fluid), public, allocatable :: tc_fl
24  !> type of fluid for thermal emission synthesis
25  type(te_fluid), public, allocatable :: te_fl_mhd
26 
27  !> Whether radiative cooling is added
28  logical, public, protected :: mhd_radiative_cooling = .false.
29  !> type of fluid for radiative cooling
30  type(rc_fluid), public, allocatable :: rc_fl
31 
32  !> Whether viscosity is added
33  logical, public, protected :: mhd_viscosity = .false.
34 
35  !> Whether gravity is added
36  logical, public, protected :: mhd_gravity = .false.
37 
38  !> Whether rotating frame is activated
39  logical, public, protected :: mhd_rotating_frame = .false.
40 
41  !> Whether Hall-MHD is used
42  logical, public, protected :: mhd_hall = .false.
43 
44  !> Whether Ambipolar term is used
45  logical, public, protected :: mhd_ambipolar = .false.
46 
47  !> Whether Ambipolar term is implemented using supertimestepping
48  logical, public, protected :: mhd_ambipolar_sts = .false.
49 
50  !> Whether Ambipolar term is implemented explicitly
51  logical, public, protected :: mhd_ambipolar_exp = .false.
52 
53  !> Whether particles module is added
54  logical, public, protected :: mhd_particles = .false.
55 
56  !> Whether magnetofriction is added
57  logical, public, protected :: mhd_magnetofriction = .false.
58 
59  !> Whether GLM-MHD is used to control div B
60  logical, public, protected :: mhd_glm = .false.
61 
62  !> Whether extended GLM-MHD is used with additional sources
63  logical, public, protected :: mhd_glm_extended = .true.
64 
65  !> Whether TRAC method is used
66  logical, public, protected :: mhd_trac = .false.
67 
68  !> Which TRAC method is used
69  integer, public, protected :: mhd_trac_type=1
70 
71  !> Height of the mask used in the TRAC method
72  double precision, public, protected :: mhd_trac_mask = 0.d0
73 
74  !> Distance between two adjacent traced magnetic field lines (in finest cell size)
75  integer, public, protected :: mhd_trac_finegrid=4
76 
77  !> Whether internal energy is solved instead of total energy
78  logical, public, protected :: mhd_internal_e = .false.
79 
80  !TODO this does not work with the splitting: check mhd_check_w_hde and mhd_handle_small_values_hde
81  !> Whether hydrodynamic energy is solved instead of total energy
82  logical, public, protected :: mhd_hydrodynamic_e = .false.
83 
84  !> Whether divB cleaning sources are added splitting from fluid solver
85  logical, public, protected :: source_split_divb = .false.
86 
87  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
88  !> taking values within [0, 1]
89  double precision, public :: mhd_glm_alpha = 0.5d0
90 
91  !TODO this does not work with the splitting: check mhd_check_w_semirelati and mhd_handle_small_values_semirelati
92  !> Whether semirelativistic MHD equations (Gombosi 2002 JCP) are solved
93  logical, public, protected :: mhd_semirelativistic = .false.
94 
95  !> Whether boris simplified semirelativistic MHD equations (Gombosi 2002 JCP) are solved
96  logical, public, protected :: mhd_boris_simplification = .false.
97 
98  !> Reduced speed of light for semirelativistic MHD: 2% of light speed
99  double precision, public, protected :: mhd_reduced_c = 0.02d0*const_c
100 
101  !> Whether plasma is partially ionized
102  logical, public, protected :: mhd_partial_ionization = .false.
103 
104  !> Whether CAK radiation line force is activated
105  logical, public, protected :: mhd_cak_force = .false.
106 
107  !> MHD fourth order
108  logical, public, protected :: mhd_4th_order = .false.
109 
110  !> whether split off equilibrium density
111  logical, public :: has_equi_rho0 = .false.
112  !> whether split off equilibrium thermal pressure
113  logical, public :: has_equi_pe0 = .false.
114  logical, public :: mhd_equi_thermal = .false.
115 
116  !> equi vars indices in the state%equi_vars array
117  integer, public :: equi_rho0_ = -1
118  integer, public :: equi_pe0_ = -1
119 
120  !> whether dump full variables (when splitting is used) in a separate dat file
121  logical, public, protected :: mhd_dump_full_vars = .false.
122 
123  !> Number of tracer species
124  integer, public, protected :: mhd_n_tracer = 0
125 
126  !> Index of the density (in the w array)
127  integer, public, protected :: rho_
128 
129  !> Indices of the momentum density
130  integer, allocatable, public, protected :: mom(:)
131 
132  !> Index of the energy density (-1 if not present)
133  integer, public, protected :: e_
134 
135  !> Index of the gas pressure (-1 if not present) should equal e_
136  integer, public, protected :: p_
137 
138 
139  !> Indices of the GLM psi
140  integer, public, protected :: psi_
141 
142  !> Indices of temperature
143  integer, public, protected :: te_
144 
145  !> Index of the cutoff temperature for the TRAC method
146  integer, public, protected :: tcoff_
147  integer, public, protected :: tweight_
148 
149  !> Indices of the tracers
150  integer, allocatable, public, protected :: tracer(:)
151 
152  !> The adiabatic index
153  double precision, public :: mhd_gamma = 5.d0/3.0d0
154 
155  !> The adiabatic constant
156  double precision, public :: mhd_adiab = 1.0d0
157 
158  !> The MHD resistivity
159  double precision, public :: mhd_eta = 0.0d0
160 
161  !> The MHD hyper-resistivity
162  double precision, public :: mhd_eta_hyper = 0.0d0
163 
164  !> TODO: what is this?
165  double precision, public :: mhd_etah = 0.0d0
166 
167  !> The MHD ambipolar coefficient
168  double precision, public :: mhd_eta_ambi = 0.0d0
169 
170  !> The small_est allowed energy
171  double precision, protected :: small_e
172 
173  !> The number of waves
174  integer :: nwwave=8
175 
176  !> Method type to clean divergence of B
177  character(len=std_len), public, protected :: typedivbfix = 'linde'
178 
179  !> Method type of constrained transport
180  character(len=std_len), public, protected :: type_ct = 'uct_contact'
181 
182  !> Whether divB is computed with a fourth order approximation
183  logical, public, protected :: mhd_divb_4thorder = .false.
184 
185  !> Method type in a integer for good performance
186  integer :: type_divb
187 
188  !> Coefficient of diffusive divB cleaning
189  double precision :: divbdiff = 0.8d0
190 
191  !> Update all equations due to divB cleaning
192  character(len=std_len) :: typedivbdiff = 'all'
193 
194  !> Use a compact way to add resistivity
195  logical :: compactres = .false.
196 
197  !> Add divB wave in Roe solver
198  logical, public :: divbwave = .true.
199 
200  !> clean initial divB
201  logical, public :: clean_initial_divb = .false.
202 
203  !> Helium abundance over Hydrogen
204  double precision, public, protected :: he_abundance=0.1d0
205  !> Ionization fraction of H
206  !> H_ion_fr = H+/(H+ + H)
207  double precision, public, protected :: h_ion_fr=1d0
208  !> Ionization fraction of He
209  !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
210  double precision, public, protected :: he_ion_fr=1d0
211  !> Ratio of number He2+ / number He+ + He2+
212  !> He_ion_fr2 = He2+/(He2+ + He+)
213  double precision, public, protected :: he_ion_fr2=1d0
214  ! used for eq of state when it is not defined by units,
215  ! the units do not contain terms related to ionization fraction
216  ! and it is p = RR * rho * T
217  double precision, public, protected :: rr=1d0
218  ! remove the below flag and assume default value = .false.
219  ! when eq state properly implemented everywhere
220  ! and not anymore through units
221  logical, public, protected :: eq_state_units = .true.
222 
223  !> To control divB=0 fix for boundary
224  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
225 
226  !> To skip * layer of ghost cells during divB=0 fix for boundary
227  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
228 
229  !> B0 field is force-free
230  logical, public, protected :: b0field_forcefree=.true.
231 
232  !> Whether an total energy equation is used
233  logical :: total_energy = .true.
234 
235  !> Whether an internal or hydrodynamic energy equation is used
236  logical :: partial_energy = .false.
237 
238  !> Whether gravity work is included in energy equation
239  logical :: gravity_energy
240 
241  !> gravity work is calculated use density times velocity or conservative momentum
242  logical :: gravity_rhov = .false.
243 
244  !> gamma minus one and its inverse
245  double precision :: gamma_1, inv_gamma_1
246 
247  !> inverse of squared speed of light c0 and reduced speed of light c
248  double precision :: inv_squared_c0, inv_squared_c
249 
250  ! DivB cleaning methods
251  integer, parameter :: divb_none = 0
252  integer, parameter :: divb_multigrid = -1
253  integer, parameter :: divb_glm = 1
254  integer, parameter :: divb_powel = 2
255  integer, parameter :: divb_janhunen = 3
256  integer, parameter :: divb_linde = 4
257  integer, parameter :: divb_lindejanhunen = 5
258  integer, parameter :: divb_lindepowel = 6
259  integer, parameter :: divb_lindeglm = 7
260  integer, parameter :: divb_ct = 8
261 
262  !define the subroutine interface for the ambipolar mask
263  abstract interface
264 
265  subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
267  integer, intent(in) :: ixi^l, ixo^l
268  double precision, intent(in) :: x(ixi^s,1:ndim)
269  double precision, intent(in) :: w(ixi^s,1:nw)
270  double precision, intent(inout) :: res(ixi^s)
271  end subroutine mask_subroutine
272 
273  function fun_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
274  use mod_global_parameters, only: nw, ndim,block
275  integer, intent(in) :: ixI^L, ixO^L
276  double precision, intent(in) :: w(ixI^S, nw)
277  double precision :: ke(ixO^S)
278  double precision, intent(in), optional :: inv_rho(ixO^S)
279  end function fun_kin_en
280 
281  end interface
282 
283  procedure(mask_subroutine), pointer :: usr_mask_ambipolar => null()
284  procedure(sub_convert), pointer :: mhd_to_primitive => null()
285  procedure(sub_convert), pointer :: mhd_to_conserved => null()
286  procedure(sub_small_values), pointer :: mhd_handle_small_values => null()
287  procedure(sub_get_pthermal), pointer :: mhd_get_pthermal => null()
288  procedure(sub_get_pthermal), pointer :: mhd_get_rfactor => null()
289  procedure(sub_get_pthermal), pointer :: mhd_get_temperature=> null()
290  procedure(sub_get_v), pointer :: mhd_get_v => null()
291  procedure(fun_kin_en), pointer :: mhd_kin_en => null()
292  ! Public methods
293  public :: usr_mask_ambipolar
294  public :: mhd_phys_init
295  public :: mhd_kin_en
296  public :: mhd_get_pthermal
297  public :: mhd_get_temperature
298  public :: mhd_get_v
299  public :: mhd_get_rho
300  public :: mhd_get_v_idim
301  public :: mhd_to_conserved
302  public :: mhd_to_primitive
303  public :: mhd_get_csound2
304  public :: mhd_e_to_ei
305  public :: mhd_ei_to_e
306  public :: mhd_face_to_center
307  public :: get_divb
308  public :: get_current
309  !> needed public if we want to use the ambipolar coefficient in the user file
310  public :: multiplyambicoef
311  public :: get_normalized_divb
312  public :: b_from_vector_potential
313  public :: mhd_mag_en_all
314  {^nooned
315  public :: mhd_clean_divb_multigrid
316  }
317 
318 contains
319 
320  !> Read this module"s parameters from a file
321  subroutine mhd_read_params(files)
323  use mod_particles, only: particles_eta, particles_etah
324  character(len=*), intent(in) :: files(:)
325  integer :: n
326 
327  namelist /mhd_list/ mhd_energy, mhd_n_tracer, mhd_gamma, mhd_adiab,&
331  typedivbdiff, type_ct, compactres, divbwave, he_abundance, &
334  particles_eta, particles_etah,has_equi_rho0, has_equi_pe0,mhd_equi_thermal,&
338 
339  do n = 1, size(files)
340  open(unitpar, file=trim(files(n)), status="old")
341  read(unitpar, mhd_list, end=111)
342 111 close(unitpar)
343  end do
344 
345  end subroutine mhd_read_params
346 
347  !> Write this module's parameters to a snapsoht
348  subroutine mhd_write_info(fh)
350  integer, intent(in) :: fh
351  integer, parameter :: n_par = 1
352  double precision :: values(n_par)
353  character(len=name_len) :: names(n_par)
354  integer, dimension(MPI_STATUS_SIZE) :: st
355  integer :: er
356 
357  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
358 
359  names(1) = "gamma"
360  values(1) = mhd_gamma
361  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
362  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
363  end subroutine mhd_write_info
364 
365  subroutine mhd_phys_init()
369  use mod_viscosity, only: viscosity_init
370  use mod_gravity, only: gravity_init
371  use mod_particles, only: particles_init, particles_eta, particles_etah
376  use mod_cak_force, only: cak_init
378  use mod_usr_methods, only: usr_rfactor
379  {^nooned
381  }
382 
383  integer :: itr, idir
384 
385  call mhd_read_params(par_files)
386 
387  if(mhd_internal_e) then
388  if(mhd_hydrodynamic_e) then
389  mhd_hydrodynamic_e=.false.
390  if(mype==0) write(*,*) 'WARNING: set mhd_hydrodynamic_e=F when mhd_internal_e=T'
391  end if
392  end if
393 
394  if(mhd_semirelativistic) then
395  if(mhd_boris_simplification) then
397  if(mype==0) write(*,*) 'WARNING: set mhd_boris_simplification=F when mhd_semirelativistic=T'
398  end if
399  if(b0field) b0fieldalloccoarse=.true.
400  end if
401 
402  if(.not. mhd_energy) then
403  if(mhd_internal_e) then
404  mhd_internal_e=.false.
405  if(mype==0) write(*,*) 'WARNING: set mhd_internal_e=F when mhd_energy=F'
406  end if
407  if(mhd_hydrodynamic_e) then
408  mhd_hydrodynamic_e=.false.
409  if(mype==0) write(*,*) 'WARNING: set mhd_hydrodynamic_e=F when mhd_energy=F'
410  end if
411  if(mhd_thermal_conduction) then
412  mhd_thermal_conduction=.false.
413  if(mype==0) write(*,*) 'WARNING: set mhd_thermal_conduction=F when mhd_energy=F'
414  end if
415  if(mhd_radiative_cooling) then
416  mhd_radiative_cooling=.false.
417  if(mype==0) write(*,*) 'WARNING: set mhd_radiative_cooling=F when mhd_energy=F'
418  end if
419  if(mhd_trac) then
420  mhd_trac=.false.
421  if(mype==0) write(*,*) 'WARNING: set mhd_trac=F when mhd_energy=F'
422  end if
423  if(mhd_partial_ionization) then
424  mhd_partial_ionization=.false.
425  if(mype==0) write(*,*) 'WARNING: set mhd_partial_ionization=F when mhd_energy=F'
426  end if
427  end if
428  if(.not.eq_state_units) then
429  if(mhd_partial_ionization) then
430  mhd_partial_ionization=.false.
431  if(mype==0) write(*,*) 'WARNING: set mhd_partial_ionization=F when eq_state_units=F'
432  end if
433  end if
434 
435  physics_type = "mhd"
436  phys_energy=mhd_energy
437  phys_internal_e=mhd_internal_e
440  phys_partial_ionization=mhd_partial_ionization
441 
442  phys_gamma = mhd_gamma
443 
444  if(mhd_energy) then
446  partial_energy=.true.
447  total_energy=.false.
448  else
449  partial_energy=.false.
450  total_energy=.true.
451  end if
452  else
453  total_energy=.false.
454  end if
455  phys_total_energy=total_energy
456  if(mhd_energy) then
457  if(mhd_internal_e) then
458  gravity_energy=.false.
459  else
460  gravity_energy=.true.
461  end if
463  gravity_rhov=.true.
464  end if
465  if(mhd_semirelativistic.and..not.mhd_hydrodynamic_e) then
466  gravity_rhov=.true.
467  end if
468  else
469  gravity_energy=.false.
470  end if
471 
472  {^ifoned
473  if(mhd_trac .and. mhd_trac_type .gt. 2) then
474  mhd_trac_type=1
475  if(mype==0) write(*,*) 'WARNING: reset mhd_trac_type=1 for 1D simulation'
476  end if
477  }
478  if(mhd_trac .and. mhd_trac_type .le. 4) then
479  mhd_trac_mask=bigdouble
480  if(mype==0) write(*,*) 'WARNING: set mhd_trac_mask==bigdouble for global TRAC method'
481  end if
483 
484  ! set default gamma for polytropic/isothermal process
486  if(ndim==1) typedivbfix='none'
487  select case (typedivbfix)
488  case ('none')
489  type_divb = divb_none
490  {^nooned
491  case ('multigrid')
492  type_divb = divb_multigrid
493  use_multigrid = .true.
494  mg%operator_type = mg_laplacian
495  phys_global_source_after => mhd_clean_divb_multigrid
496  }
497  case ('glm')
498  mhd_glm = .true.
499  need_global_cmax = .true.
500  type_divb = divb_glm
501  case ('powel', 'powell')
502  type_divb = divb_powel
503  case ('janhunen')
504  type_divb = divb_janhunen
505  case ('linde')
506  type_divb = divb_linde
507  case ('lindejanhunen')
508  type_divb = divb_lindejanhunen
509  case ('lindepowel')
510  type_divb = divb_lindepowel
511  case ('lindeglm')
512  mhd_glm = .true.
513  need_global_cmax = .true.
514  type_divb = divb_lindeglm
515  case ('ct')
516  type_divb = divb_ct
517  stagger_grid = .true.
518  case default
519  call mpistop('Unknown divB fix')
520  end select
521 
522  allocate(start_indices(number_species),stop_indices(number_species))
523  ! set the index of the first flux variable for species 1
524  start_indices(1)=1
525  ! Determine flux variables
526  rho_ = var_set_rho()
527 
528  allocate(mom(ndir))
529  mom(:) = var_set_momentum(ndir)
530 
531  ! Set index of energy variable
532  if (mhd_energy) then
533  nwwave = 8
534  e_ = var_set_energy() ! energy density
535  p_ = e_ ! gas pressure
536  else
537  nwwave = 7
538  e_ = -1
539  p_ = -1
540  end if
541 
542  allocate(mag(ndir))
543  mag(:) = var_set_bfield(ndir)
544 
545  if (mhd_glm) then
546  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
547  else
548  psi_ = -1
549  end if
550 
551  allocate(tracer(mhd_n_tracer))
552  ! Set starting index of tracers
553  do itr = 1, mhd_n_tracer
554  tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
555  end do
556 
557  ! set number of variables which need update ghostcells
558  nwgc=nwflux
559 
560  ! set the index of the last flux variable for species 1
561  stop_indices(1)=nwflux
562 
563  ! set temperature as an auxiliary variable to get ionization degree
564  if(mhd_partial_ionization) then
565  te_ = var_set_auxvar('Te','Te')
566  else
567  te_ = -1
568  end if
569 
570  ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
571  tweight_ = -1
572  if(mhd_trac) then
573  tcoff_ = var_set_wextra()
574  iw_tcoff=tcoff_
575  if(mhd_trac_type .ge. 3) then
576  tweight_ = var_set_wextra()
577  endif
578  else
579  tcoff_ = -1
580  end if
581 
582  ! set indices of equi vars and update number_equi_vars
583  number_equi_vars = 0
584  if(has_equi_rho0) then
587  iw_equi_rho = equi_rho0_
588  endif
589  if(has_equi_pe0) then
592  iw_equi_p = equi_pe0_
593  phys_equi_pe=.true.
594  endif
595  ! determine number of stagger variables
596  nws=ndim
597 
598  nvector = 2 ! No. vector vars
599  allocate(iw_vector(nvector))
600  iw_vector(1) = mom(1) - 1 ! TODO: why like this?
601  iw_vector(2) = mag(1) - 1 ! TODO: why like this?
602 
603  ! Check whether custom flux types have been defined
604  if (.not. allocated(flux_type)) then
605  allocate(flux_type(ndir, nwflux))
606  flux_type = flux_default
607  else if (any(shape(flux_type) /= [ndir, nwflux])) then
608  call mpistop("phys_check error: flux_type has wrong shape")
609  end if
610 
611  if(nwflux>mag(ndir)) then
612  ! for flux of tracers, using hll flux
613  flux_type(:,mag(ndir)+1:nwflux)=flux_hll
614  end if
615 
616  if(ndim>1) then
617  if(mhd_glm) then
618  flux_type(:,psi_)=flux_special
619  do idir=1,ndir
620  flux_type(idir,mag(idir))=flux_special
621  end do
622  else
623  do idir=1,ndir
624  flux_type(idir,mag(idir))=flux_tvdlf
625  end do
626  end if
627  end if
628 
629  phys_get_dt => mhd_get_dt
630  if(mhd_semirelativistic) then
631  phys_get_cmax => mhd_get_cmax_semirelati
632  else
633  phys_get_cmax => mhd_get_cmax_origin
634  end if
635  phys_get_a2max => mhd_get_a2max
636  phys_get_tcutoff => mhd_get_tcutoff
637  phys_get_h_speed => mhd_get_h_speed
638  if(has_equi_rho0) then
639  phys_get_cbounds => mhd_get_cbounds_split_rho
640  else if(mhd_semirelativistic) then
641  phys_get_cbounds => mhd_get_cbounds_semirelati
642  else
643  phys_get_cbounds => mhd_get_cbounds
644  end if
645  if(mhd_hydrodynamic_e) then
646  phys_to_primitive => mhd_to_primitive_hde
648  phys_to_conserved => mhd_to_conserved_hde
650  else if(mhd_semirelativistic) then
651  phys_to_primitive => mhd_to_primitive_semirelati
653  phys_to_conserved => mhd_to_conserved_semirelati
655  else
656  if(has_equi_rho0) then
657  phys_to_primitive => mhd_to_primitive_split_rho
659  phys_to_conserved => mhd_to_conserved_split_rho
661  else if(mhd_internal_e) then
662  phys_to_primitive => mhd_to_primitive_inte
664  phys_to_conserved => mhd_to_conserved_inte
666  else
667  phys_to_primitive => mhd_to_primitive_origin
669  phys_to_conserved => mhd_to_conserved_origin
671  end if
672  end if
673  if(mhd_hydrodynamic_e) then
674  phys_get_flux => mhd_get_flux_hde
675  else if(mhd_semirelativistic) then
676  phys_get_flux => mhd_get_flux_semirelati
677  else
678  if(b0field.or.has_equi_rho0.or.has_equi_pe0) then
679  phys_get_flux => mhd_get_flux_split
680  else
681  phys_get_flux => mhd_get_flux
682  end if
683  end if
684  if(mhd_boris_simplification) then
685  phys_get_v => mhd_get_v_boris
688  else
689  phys_get_v => mhd_get_v_origin
692  end if
693  if(b0field.or.has_equi_rho0) then
694  phys_add_source_geom => mhd_add_source_geom_split
695  else
696  phys_add_source_geom => mhd_add_source_geom
697  end if
698  phys_add_source => mhd_add_source
699  phys_check_params => mhd_check_params
700  phys_write_info => mhd_write_info
701 
702  if(mhd_hydrodynamic_e) then
703  phys_handle_small_values => mhd_handle_small_values_hde
704  mhd_handle_small_values => mhd_handle_small_values_hde
705  phys_check_w => mhd_check_w_hde
706  else if(mhd_semirelativistic) then
707  phys_handle_small_values => mhd_handle_small_values_semirelati
708  mhd_handle_small_values => mhd_handle_small_values_semirelati
709  phys_check_w => mhd_check_w_semirelati
710  else
711  if(mhd_internal_e) then
712  phys_handle_small_values => mhd_handle_small_values_inte
713  mhd_handle_small_values => mhd_handle_small_values_inte
714  phys_check_w => mhd_check_w_inte
715  else
716  phys_handle_small_values => mhd_handle_small_values_origin
717  mhd_handle_small_values => mhd_handle_small_values_origin
718  phys_check_w => mhd_check_w_origin
719  end if
720  end if
721 
722  if(.not.mhd_energy) then
723  phys_get_pthermal => mhd_get_pthermal_iso
725  else
726  if(mhd_internal_e) then
727  phys_get_pthermal => mhd_get_pthermal_eint
729  else if(mhd_hydrodynamic_e) then
730  phys_get_pthermal => mhd_get_pthermal_hde
732  else if(mhd_semirelativistic) then
733  phys_get_pthermal => mhd_get_pthermal_semirelati
735  else
736  phys_get_pthermal => mhd_get_pthermal_origin
738  end if
739  end if
740  if(number_equi_vars>0) then
741  phys_set_equi_vars => set_equi_vars_grid
742  endif
743 
744  if(type_divb==divb_glm) then
745  phys_modify_wlr => mhd_modify_wlr
746  end if
747 
748  ! choose Rfactor in ideal gas law
749  if(mhd_partial_ionization) then
750  mhd_get_rfactor=>rfactor_from_temperature_ionization
751  phys_update_temperature => mhd_update_temperature
752  else if(associated(usr_rfactor)) then
753  mhd_get_rfactor=>usr_rfactor
754  else
755  mhd_get_rfactor=>rfactor_from_constant_ionization
756  end if
757 
758  if(mhd_partial_ionization) then
760  else
761  if(mhd_internal_e) then
762  if(has_equi_pe0 .and. has_equi_rho0) then
764  else
766  end if
767  else
768  if(has_equi_pe0 .and. has_equi_rho0) then
770  else
772  end if
773  end if
774  end if
775 
776  ! if using ct stagger grid, boundary divb=0 is not done here
777  if(stagger_grid) then
778  phys_get_ct_velocity => mhd_get_ct_velocity
779  phys_update_faces => mhd_update_faces
780  phys_face_to_center => mhd_face_to_center
781  phys_modify_wlr => mhd_modify_wlr
782  else if(ndim>1) then
783  phys_boundary_adjust => mhd_boundary_adjust
784  end if
785 
786  {^nooned
787  ! clean initial divb
788  if(clean_initial_divb) phys_clean_divb => mhd_clean_divb_multigrid
789  }
790 
791  ! Whether diagonal ghost cells are required for the physics
792  if(type_divb < divb_linde) phys_req_diagonal = .false.
793 
794  ! derive units from basic units
795  call mhd_physical_units()
796 
797  if(.not. mhd_energy .and. mhd_thermal_conduction) then
798  call mpistop("thermal conduction needs mhd_energy=T")
799  end if
800  if(.not. mhd_energy .and. mhd_radiative_cooling) then
801  call mpistop("radiative cooling needs mhd_energy=T")
802  end if
803 
804  ! resistive MHD needs diagonal ghost cells
805  if(mhd_eta/=0.d0) phys_req_diagonal = .true.
806 
807  ! initialize thermal conduction module
808  if (mhd_thermal_conduction) then
809  phys_req_diagonal = .true.
810 
811  call sts_init()
813 
814  allocate(tc_fl)
817  if(phys_internal_e) then
818  if(has_equi_pe0 .and. has_equi_rho0) then
819  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint_with_equi
820  else
821  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_eint
822  end if
823  else
824  if(has_equi_pe0 .and. has_equi_rho0) then
825  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot_with_equi
826  else
827  tc_fl%get_temperature_from_conserved => mhd_get_temperature_from_etot
828  end if
829  end if
830  if(has_equi_pe0 .and. has_equi_rho0) then
831  tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
832  if(mhd_equi_thermal) then
833  tc_fl%has_equi = .true.
834  tc_fl%get_temperature_equi => mhd_get_temperature_equi
835  tc_fl%get_rho_equi => mhd_get_rho_equi
836  else
837  tc_fl%has_equi = .false.
838  end if
839  else
840  tc_fl%get_temperature_from_eint => mhd_get_temperature_from_eint
841  end if
842  if(.not.mhd_internal_e) then
843  if(mhd_hydrodynamic_e) then
845  else if(mhd_semirelativistic) then
847  else
849  end if
850  end if
852  tc_fl%get_rho => mhd_get_rho
853  tc_fl%e_ = e_
854  tc_fl%Tcoff_ = tcoff_
855  end if
856 
857  ! Initialize radiative cooling module
858  if (mhd_radiative_cooling) then
860  allocate(rc_fl)
862  rc_fl%get_rho => mhd_get_rho
863  rc_fl%get_pthermal => mhd_get_pthermal
864  rc_fl%get_var_Rfactor => mhd_get_rfactor
865  rc_fl%e_ = e_
866  rc_fl%Tcoff_ = tcoff_
867  if(has_equi_pe0 .and. has_equi_rho0 .and. mhd_equi_thermal) then
868  rc_fl%has_equi = .true.
869  rc_fl%get_rho_equi => mhd_get_rho_equi
870  rc_fl%get_pthermal_equi => mhd_get_pe_equi
871  else
872  rc_fl%has_equi = .false.
873  end if
874  end if
875  allocate(te_fl_mhd)
876  te_fl_mhd%get_rho=> mhd_get_rho
877  te_fl_mhd%get_pthermal=> mhd_get_pthermal
878  te_fl_mhd%get_var_Rfactor => mhd_get_rfactor
879 {^ifthreed
880  phys_te_images => mhd_te_images
881 }
882  ! Initialize viscosity module
883  if (mhd_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
884 
885  ! Initialize gravity module
886  if(mhd_gravity) then
887  call gravity_init()
888  end if
889 
890  ! Initialize rotating frame module
892 
893  ! Initialize particles module
894  if(mhd_particles) then
895  call particles_init()
896  if (particles_eta < zero) particles_eta = mhd_eta
897  if (particles_etah < zero) particles_eta = mhd_etah
898  phys_req_diagonal = .true.
899  if(mype==0) then
900  write(*,*) '*****Using particles: with mhd_eta, mhd_etah :', mhd_eta, mhd_etah
901  write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
902  end if
903  end if
904 
905  ! initialize magnetofriction module
906  if(mhd_magnetofriction) then
907  phys_req_diagonal = .true.
908  call magnetofriction_init()
909  end if
910 
911  ! For Hall, we need one more reconstructed layer since currents are computed
912  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
913  ! added in nghostcells.
914  if(mhd_hall) then
915  phys_req_diagonal = .true.
916  if(mhd_4th_order) then
917  phys_wider_stencil = 2
918  else
919  phys_wider_stencil = 1
920  end if
921  end if
922 
923  if(mhd_ambipolar) then
924  phys_req_diagonal = .true.
925  if(mhd_ambipolar_sts) then
926  call sts_init()
927  if(mhd_internal_e) then
929  ndir,mag(1),ndir,.true.)
930  else
932  mag(ndir)-mom(ndir),mag(1),ndir,.true.)
933  end if
934  else
935  mhd_ambipolar_exp=.true.
936  ! For flux ambipolar term, we need one more reconstructed layer since currents are computed
937  ! in mhd_get_flux: assuming one additional ghost layer (two for FOURTHORDER) was
938  ! added in nghostcells.
939  if(mhd_4th_order) then
940  phys_wider_stencil = 2
941  else
942  phys_wider_stencil = 1
943  end if
944  end if
945  end if
946 
947  ! initialize ionization degree table
949 
950  ! Initialize CAK radiation force module
951  if (mhd_cak_force) call cak_init(mhd_gamma)
952 
953  end subroutine mhd_phys_init
954 
955 {^ifthreed
956  subroutine mhd_te_images
959 
960  select case(convert_type)
961  case('EIvtiCCmpi','EIvtuCCmpi')
963  case('ESvtiCCmpi','ESvtuCCmpi')
965  case('SIvtiCCmpi','SIvtuCCmpi')
967  case('WIvtiCCmpi','WIvtuCCmpi')
969  case default
970  call mpistop("Error in synthesize emission: Unknown convert_type")
971  end select
972  end subroutine mhd_te_images
973 }
974 
975 !!start th cond
976  ! wrappers for STS functions in thermal_conductivity module
977  ! which take as argument the tc_fluid (defined in the physics module)
978  subroutine mhd_sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
980  use mod_fix_conserve
982  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
983  double precision, intent(in) :: x(ixI^S,1:ndim)
984  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
985  double precision, intent(in) :: my_dt
986  logical, intent(in) :: fix_conserve_at_step
987  call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
988  end subroutine mhd_sts_set_source_tc_mhd
989 
990  function mhd_get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
991  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
992  !where tc_k_para_i=tc_k_para*B_i**2/B**2
993  !and T=p/rho
996 
997  integer, intent(in) :: ixi^l, ixo^l
998  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
999  double precision, intent(in) :: w(ixi^s,1:nw)
1000  double precision :: dtnew
1001 
1002  dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
1003  end function mhd_get_tc_dt_mhd
1004 
1005  subroutine mhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
1007 
1008  integer, intent(in) :: ixI^L,ixO^L
1009  double precision, intent(inout) :: w(ixI^S,1:nw)
1010  double precision, intent(in) :: x(ixI^S,1:ndim)
1011  integer, intent(in) :: step
1012  character(len=140) :: error_msg
1013 
1014  write(error_msg,"(a,i3)") "Thermal conduction step ", step
1015  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
1016  end subroutine mhd_tc_handle_small_e
1017 
1018  ! fill in tc_fluid fields from namelist
1019  subroutine tc_params_read_mhd(fl)
1021  type(tc_fluid), intent(inout) :: fl
1022 
1023  integer :: n
1024 
1025  ! list parameters
1026  logical :: tc_perpendicular=.false.
1027  logical :: tc_saturate=.false.
1028  double precision :: tc_k_para=0d0
1029  double precision :: tc_k_perp=0d0
1030  character(len=std_len) :: tc_slope_limiter="MC"
1031 
1032  namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1033 
1034  do n = 1, size(par_files)
1035  open(unitpar, file=trim(par_files(n)), status="old")
1036  read(unitpar, tc_list, end=111)
1037 111 close(unitpar)
1038  end do
1039 
1040  fl%tc_perpendicular = tc_perpendicular
1041  fl%tc_saturate = tc_saturate
1042  fl%tc_k_para = tc_k_para
1043  fl%tc_k_perp = tc_k_perp
1044  select case(tc_slope_limiter)
1045  case ('no','none')
1046  fl%tc_slope_limiter = 0
1047  case ('MC')
1048  ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
1049  fl%tc_slope_limiter = 1
1050  case('minmod')
1051  ! minmod limiter
1052  fl%tc_slope_limiter = 2
1053  case ('superbee')
1054  ! Roes superbee limiter (eq.3.51i)
1055  fl%tc_slope_limiter = 3
1056  case ('koren')
1057  ! Barry Koren Right variant
1058  fl%tc_slope_limiter = 4
1059  case default
1060  call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
1061  end select
1062  end subroutine tc_params_read_mhd
1063 !!end th cond
1064 
1065 !!rad cool
1066  subroutine rc_params_read(fl)
1068  use mod_constants, only: bigdouble
1069  type(rc_fluid), intent(inout) :: fl
1070  integer :: n
1071  ! list parameters
1072  integer :: ncool = 4000
1073  double precision :: cfrac=0.1d0
1074 
1075  !> Name of cooling curve
1076  character(len=std_len) :: coolcurve='JCcorona'
1077 
1078  !> Name of cooling method
1079  character(len=std_len) :: coolmethod='exact'
1080 
1081  !> Fixed temperature not lower than tlow
1082  logical :: Tfix=.false.
1083 
1084  !> Lower limit of temperature
1085  double precision :: tlow=bigdouble
1086 
1087  !> Add cooling source in a split way (.true.) or un-split way (.false.)
1088  logical :: rc_split=.false.
1089 
1090 
1091  namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1092 
1093  do n = 1, size(par_files)
1094  open(unitpar, file=trim(par_files(n)), status="old")
1095  read(unitpar, rc_list, end=111)
1096 111 close(unitpar)
1097  end do
1098 
1099  fl%ncool=ncool
1100  fl%coolcurve=coolcurve
1101  fl%coolmethod=coolmethod
1102  fl%tlow=tlow
1103  fl%Tfix=tfix
1104  fl%rc_split=rc_split
1105  fl%cfrac=cfrac
1106 
1107  end subroutine rc_params_read
1108 !! end rad cool
1109 
1110  !> sets the equilibrium variables
1111  subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1113  use mod_usr_methods
1114  integer, intent(in) :: igrid, ixI^L, ixO^L
1115  double precision, intent(in) :: x(ixI^S,1:ndim)
1116 
1117  double precision :: delx(ixI^S,1:ndim)
1118  double precision :: xC(ixI^S,1:ndim),xshift^D
1119  integer :: idims, ixC^L, hxO^L, ix, idims2
1120 
1121  if(slab_uniform)then
1122  ^d&delx(ixi^s,^d)=rnode(rpdx^d_,igrid)\
1123  else
1124  ! for all non-cartesian and stretched cartesian coordinates
1125  delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1126  endif
1127 
1128  do idims=1,ndim
1129  hxo^l=ixo^l-kr(idims,^d);
1130  if(stagger_grid) then
1131  ! ct needs all transverse cells
1132  ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
1133  else
1134  ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
1135  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1136  end if
1137  ! always xshift=0 or 1/2
1138  xshift^d=half*(one-kr(^d,idims));
1139  do idims2=1,ndim
1140  select case(idims2)
1141  {case(^d)
1142  do ix = ixc^lim^d
1143  ! xshift=half: this is the cell center coordinate
1144  ! xshift=0: this is the cell edge i+1/2 coordinate
1145  xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1146  end do\}
1147  end select
1148  end do
1149  call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1150  end do
1151 
1152  end subroutine set_equi_vars_grid_faces
1153 
1154  !> sets the equilibrium variables
1155  subroutine set_equi_vars_grid(igrid)
1157  use mod_usr_methods
1158 
1159  integer, intent(in) :: igrid
1160 
1161  !values at the center
1162  call usr_set_equi_vars(ixg^ll,ixg^ll,ps(igrid)%x,ps(igrid)%equi_vars(ixg^t,1:number_equi_vars,0))
1163 
1164  !values at the interfaces
1165  call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^ll,ixm^ll)
1166 
1167  end subroutine set_equi_vars_grid
1168 
1169  ! w, wnew conserved, add splitted variables back to wnew
1170  function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc) result(wnew)
1172  integer, intent(in) :: ixi^l,ixo^l, nwc
1173  double precision, intent(in) :: w(ixi^s, 1:nw)
1174  double precision, intent(in) :: x(ixi^s,1:ndim)
1175  double precision :: wnew(ixo^s, 1:nwc)
1176  double precision :: rho(ixi^s)
1177 
1178  call mhd_get_rho(w,x,ixi^l,ixo^l,rho(ixi^s))
1179  wnew(ixo^s,rho_) = rho(ixo^s)
1180  wnew(ixo^s,mom(:)) = w(ixo^s,mom(:))
1181 
1182  if (b0field) then
1183  ! add background magnetic field B0 to B
1184  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+block%B0(ixo^s,:,0)
1185  else
1186  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1187  end if
1188 
1189  if(mhd_energy) then
1190  wnew(ixo^s,e_) = w(ixo^s,e_)
1191  if(has_equi_pe0) then
1192  wnew(ixo^s,e_) = wnew(ixo^s,e_) + block%equi_vars(ixo^s,equi_pe0_,0)* inv_gamma_1
1193  end if
1194  if(b0field .and. .not. mhd_internal_e) then
1195  wnew(ixo^s,e_)=wnew(ixo^s,e_)+0.5d0*sum(block%B0(ixo^s,:,0)**2,dim=ndim+1) &
1196  + sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,0),dim=ndim+1)
1197  end if
1198  end if
1199 
1200  end function convert_vars_splitting
1201 
1202  subroutine mhd_check_params
1204  use mod_usr_methods
1205  use mod_convert, only: add_convert_method
1206 
1207  ! after user parameter setting
1208  gamma_1=mhd_gamma-1.d0
1209  if (.not. mhd_energy) then
1210  if (mhd_gamma <= 0.0d0) call mpistop ("Error: mhd_gamma <= 0")
1211  if (mhd_adiab < 0.0d0) call mpistop ("Error: mhd_adiab < 0")
1213  else
1214  if (mhd_gamma <= 0.0d0 .or. mhd_gamma == 1.0d0) &
1215  call mpistop ("Error: mhd_gamma <= 0 or mhd_gamma == 1")
1216  inv_gamma_1=1.d0/gamma_1
1217  small_e = small_pressure * inv_gamma_1
1218  end if
1219 
1220  if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
1221  call mpistop("usr_set_equi_vars has to be implemented in the user file")
1222  endif
1223  if(convert .or. autoconvert) then
1224  if(convert_type .eq. 'dat_generic_mpi') then
1225  if(mhd_dump_full_vars) then
1226  if(mype .eq. 0) print*, " add conversion method: split -> full "
1227  call add_convert_method(convert_vars_splitting, nw, cons_wnames, "new")
1228  endif
1229  endif
1230  endif
1231  end subroutine mhd_check_params
1232 
1233  subroutine mhd_physical_units()
1235  double precision :: mp,kB,miu0,c_lightspeed
1236  double precision :: a,b
1237  ! Derive scaling units
1238  if(si_unit) then
1239  mp=mp_si
1240  kb=kb_si
1241  miu0=miu0_si
1242  c_lightspeed=c_si
1243  else
1244  mp=mp_cgs
1245  kb=kb_cgs
1246  miu0=4.d0*dpi ! G^2 cm^2 dyne^-1
1247  c_lightspeed=const_c
1248  end if
1249  if(eq_state_units) then
1250  a = 1d0 + 4d0 * he_abundance
1251  if(mhd_partial_ionization) then
1252  b = 2+.3d0
1253  else
1254  b = 1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0)
1255  end if
1256  rr = 1d0
1257  else
1258  a = 1d0
1259  b = 1d0
1260  rr = (1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0))/(1d0 + 4d0 * he_abundance)
1261  end if
1262  if(unit_density/=1.d0) then
1264  else
1265  ! unit of numberdensity is independent by default
1267  end if
1268  if(unit_velocity/=1.d0) then
1272  else if(unit_pressure/=1.d0) then
1276  else if(unit_magneticfield/=1.d0) then
1280  else if(unit_temperature/=1.d0) then
1284  end if
1285  if(unit_time/=1.d0) then
1287  else
1288  ! unit of length is independent by default
1290  end if
1291  ! Additional units needed for the particles
1292  c_norm=c_lightspeed/unit_velocity
1294  if (.not. si_unit) unit_charge = unit_charge*const_c
1296 
1298  if(mhd_reduced_c<1.d0) then
1299  ! dimensionless speed
1300  inv_squared_c0=1.d0
1301  inv_squared_c=1.d0/mhd_reduced_c**2
1302  else
1303  inv_squared_c0=(unit_velocity/c_lightspeed)**2
1304  inv_squared_c=(unit_velocity/mhd_reduced_c)**2
1305  end if
1306  end if
1307 
1308  end subroutine mhd_physical_units
1309 
1310  subroutine mhd_check_w_semirelati(primitive,ixI^L,ixO^L,w,flag)
1312 
1313  logical, intent(in) :: primitive
1314  integer, intent(in) :: ixI^L, ixO^L
1315  double precision, intent(in) :: w(ixI^S,nw)
1316  double precision :: tmp(ixO^S),b2(ixO^S),b(ixO^S,1:ndir),Ba(ixO^S,1:ndir)
1317  double precision :: v(ixO^S,1:ndir),gamma2(ixO^S),inv_rho(ixO^S)
1318  logical, intent(inout) :: flag(ixI^S,1:nw)
1319 
1320  integer :: idir, jdir, kdir
1321 
1322  flag=.false.
1323  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1324 
1325  if(mhd_energy) then
1326  if(primitive) then
1327  where(w(ixo^s,p_) < small_pressure) flag(ixo^s,e_) = .true.
1328  else
1329  if(mhd_internal_e) then
1330  tmp(ixo^s)=w(ixo^s,e_)
1331  else
1332  if(b0field) then
1333  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1334  else
1335  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1336  end if
1337  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1338  b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1339  tmp(ixo^s)=sqrt(b2(ixo^s))
1340  where(tmp(ixo^s)>smalldouble)
1341  tmp(ixo^s)=1.d0/tmp(ixo^s)
1342  else where
1343  tmp(ixo^s)=0.d0
1344  end where
1345  do idir=1,ndir
1346  b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1347  end do
1348  tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,mom(:)),dim=ndim+1)
1349  ! Va^2/c^2
1350  b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1351  ! equation (15)
1352  gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1353  ! Convert momentum to velocity
1354  do idir = 1, ndir
1355  v(ixo^s,idir) = gamma2*(w(ixo^s, mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1356  end do
1357  ! E=Bxv
1358  b=0.d0
1359  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1360  if(lvc(idir,jdir,kdir)==1)then
1361  b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*v(ixo^s,kdir)
1362  else if(lvc(idir,jdir,kdir)==-1)then
1363  b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*v(ixo^s,kdir)
1364  end if
1365  end do; end do; end do
1366  ! Calculate internal e = e-eK-eB-eE
1367  tmp(ixo^s)=w(ixo^s,e_)&
1368  -half*(sum(v(ixo^s,:)**2,dim=ndim+1)*w(ixo^s,rho_)&
1369  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1370  +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1371  end if
1372  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1373  end if
1374  end if
1375 
1376  end subroutine mhd_check_w_semirelati
1377 
1378  subroutine mhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
1380 
1381  logical, intent(in) :: primitive
1382  integer, intent(in) :: ixI^L, ixO^L
1383  double precision, intent(in) :: w(ixI^S,nw)
1384  double precision :: tmp(ixI^S)
1385  logical, intent(inout) :: flag(ixI^S,1:nw)
1386 
1387  flag=.false.
1388  if(has_equi_rho0) then
1389  tmp(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
1390  where(tmp(ixo^s) < small_density) flag(ixo^s,rho_) = .true.
1391  else
1392  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1393  end if
1394 
1395  if(mhd_energy) then
1396  if(primitive) then
1397  if(has_equi_pe0) then
1398  tmp(ixo^s) = w(ixo^s,e_)+block%equi_vars(ixo^s,equi_pe0_,0)
1399  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
1400  else
1401  where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
1402  end if
1403  else
1404  tmp(ixo^s)=w(ixo^s,e_)-&
1405  mhd_kin_en(w,ixi^l,ixo^l)-mhd_mag_en(w,ixi^l,ixo^l)
1406  if(has_equi_pe0) then
1407  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
1408  end if
1409  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1410  end if
1411  end if
1412 
1413  end subroutine mhd_check_w_origin
1414 
1415  subroutine mhd_check_w_inte(primitive,ixI^L,ixO^L,w,flag)
1417 
1418  logical, intent(in) :: primitive
1419  integer, intent(in) :: ixI^L, ixO^L
1420  double precision, intent(in) :: w(ixI^S,nw)
1421  double precision :: tmp(ixI^S)
1422  logical, intent(inout) :: flag(ixI^S,1:nw)
1423 
1424  flag=.false.
1425  if(has_equi_rho0) then
1426  tmp(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
1427  where(tmp(ixo^s) < small_density) flag(ixo^s,rho_) = .true.
1428  else
1429  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1430  end if
1431 
1432  if(mhd_energy) then
1433  if(primitive) then
1434  if(has_equi_pe0) then
1435  tmp(ixo^s) = w(ixo^s,e_)+block%equi_vars(ixo^s,equi_pe0_,0)
1436  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_) = .true.
1437  else
1438  where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
1439  end if
1440  else
1441  if(has_equi_pe0) then
1442  tmp(ixo^s) = w(ixo^s,e_)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
1443  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1444  else
1445  where(w(ixo^s,e_) < small_e) flag(ixo^s,e_) = .true.
1446  end if
1447  end if
1448  end if
1449 
1450  end subroutine mhd_check_w_inte
1451 
1452  subroutine mhd_check_w_hde(primitive,ixI^L,ixO^L,w,flag)
1454 
1455  logical, intent(in) :: primitive
1456  integer, intent(in) :: ixI^L, ixO^L
1457  double precision, intent(in) :: w(ixI^S,nw)
1458  double precision :: tmp(ixI^S)
1459  logical, intent(inout) :: flag(ixI^S,1:nw)
1460 
1461  flag=.false.
1462  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1463 
1464  if(mhd_energy) then
1465  if(primitive) then
1466  where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
1467  else
1468  tmp(ixo^s)=w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l)
1469  where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
1470  end if
1471  end if
1472 
1473  end subroutine mhd_check_w_hde
1474 
1475  !> Transform primitive variables into conservative ones
1476  subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
1478  integer, intent(in) :: ixI^L, ixO^L
1479  double precision, intent(inout) :: w(ixI^S, nw)
1480  double precision, intent(in) :: x(ixI^S, 1:ndim)
1481 
1482  double precision :: inv_gamma2(ixO^S)
1483  integer :: idir
1484 
1485  !if (fix_small_values) then
1486  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1487  !end if
1488 
1489  ! Calculate total energy from pressure, kinetic and magnetic energy
1490  if(mhd_energy) then
1491  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1492  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1493  +mhd_mag_en(w, ixi^l, ixo^l)
1494  end if
1495 
1496  if(mhd_boris_simplification) then
1497  ! rho*(1+B^2/rho/c^2)
1498  inv_gamma2=w(ixo^s,rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1499  ! Convert velocity to momentum
1500  do idir = 1, ndir
1501  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s, mom(idir))
1502  end do
1503  else
1504  ! Convert velocity to momentum
1505  do idir = 1, ndir
1506  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1507  end do
1508  end if
1509  end subroutine mhd_to_conserved_origin
1510 
1511  !> Transform primitive variables into conservative ones
1512  subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
1514  integer, intent(in) :: ixI^L, ixO^L
1515  double precision, intent(inout) :: w(ixI^S, nw)
1516  double precision, intent(in) :: x(ixI^S, 1:ndim)
1517 
1518  double precision :: inv_gamma2(ixO^S)
1519  integer :: idir
1520 
1521  ! Calculate total energy from pressure, kinetic and magnetic energy
1522  if(mhd_energy) then
1523  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1524  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)
1525  end if
1526 
1527  ! Convert velocity to momentum
1528  if(mhd_boris_simplification) then
1529  ! rho*(1+B^2/rho/c^2)
1530  inv_gamma2=w(ixo^s,rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1531  ! Convert velocity to momentum
1532  do idir = 1, ndir
1533  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s, mom(idir))
1534  end do
1535  else
1536  do idir = 1, ndir
1537  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1538  end do
1539  end if
1540  end subroutine mhd_to_conserved_hde
1541 
1542  !> Transform primitive variables into conservative ones
1543  subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
1545  integer, intent(in) :: ixI^L, ixO^L
1546  double precision, intent(inout) :: w(ixI^S, nw)
1547  double precision, intent(in) :: x(ixI^S, 1:ndim)
1548 
1549  double precision :: inv_gamma2(ixO^S)
1550  integer :: idir
1551 
1552  ! Calculate total energy from pressure, kinetic and magnetic energy
1553  if(mhd_energy) then
1554  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
1555  end if
1556 
1557  if(mhd_boris_simplification) then
1558  ! rho*(1+B^2/rho/c^2)
1559  inv_gamma2=w(ixo^s,rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1560  ! Convert velocity to momentum
1561  do idir = 1, ndir
1562  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s, mom(idir))
1563  end do
1564  else
1565  ! Convert velocity to momentum
1566  do idir = 1, ndir
1567  w(ixo^s, mom(idir)) = w(ixo^s,rho_)*w(ixo^s, mom(idir))
1568  end do
1569  end if
1570  end subroutine mhd_to_conserved_inte
1571 
1572  !> Transform primitive variables into conservative ones
1573  subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
1575  integer, intent(in) :: ixI^L, ixO^L
1576  double precision, intent(inout) :: w(ixI^S, nw)
1577  double precision, intent(in) :: x(ixI^S, 1:ndim)
1578 
1579  double precision :: rho(ixI^S), inv_gamma2(ixO^S)
1580  integer :: idir
1581 
1582  !if (fix_small_values) then
1583  ! call mhd_handle_small_values(.true., w, x, ixI^L, ixO^L, 'mhd_to_conserved')
1584  !end if
1585 
1586  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
1587  ! Calculate total energy from pressure, kinetic and magnetic energy
1588  if(mhd_energy) then
1589  if(mhd_internal_e) then
1590  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
1591  else
1592  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1593  +half*sum(w(ixo^s,mom(:))**2,dim=ndim+1)*rho(ixo^s)&
1594  +mhd_mag_en(w, ixi^l, ixo^l)
1595  end if
1596  end if
1597 
1598  if(mhd_boris_simplification) then
1599  ! rho*(1+B^2/rho/c^2)
1600  inv_gamma2=w(ixo^s,rho_)+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_squared_c
1601  ! Convert velocity to momentum
1602  do idir = 1, ndir
1603  w(ixo^s, mom(idir)) = inv_gamma2*w(ixo^s, mom(idir))
1604  end do
1605  else
1606  ! Convert velocity to momentum
1607  do idir = 1, ndir
1608  w(ixo^s, mom(idir)) = rho(ixo^s) * w(ixo^s, mom(idir))
1609  end do
1610  end if
1611  end subroutine mhd_to_conserved_split_rho
1612 
1613  !> Transform primitive variables into conservative ones
1614  subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
1616  integer, intent(in) :: ixI^L, ixO^L
1617  double precision, intent(inout) :: w(ixI^S, nw)
1618  double precision, intent(in) :: x(ixI^S, 1:ndim)
1619 
1620  double precision :: E(ixO^S,1:ndir), B(ixO^S,1:ndir), S(ixO^S,1:ndir)
1621  integer :: idir, jdir, kdir
1622 
1623  if(b0field) then
1624  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1625  else
1626  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1627  end if
1628  e=0.d0
1629  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1630  if(lvc(idir,jdir,kdir)==1)then
1631  e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,mom(kdir))
1632  else if(lvc(idir,jdir,kdir)==-1)then
1633  e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,mom(kdir))
1634  end if
1635  end do; end do; end do
1636 
1637  if(mhd_internal_e) then
1638  ! internal energy
1639  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1
1640  else if(mhd_energy) then
1641  ! equation (9)
1642  ! Calculate total energy from internal, kinetic and magnetic energy
1643  w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1&
1644  +half*(sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1645  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1646  +sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
1647  end if
1648 
1649  ! Convert velocity to momentum
1650  do idir = 1, ndir
1651  w(ixo^s, mom(idir)) = w(ixo^s,rho_) * w(ixo^s, mom(idir))
1652  end do
1653  ! equation (5) Poynting vector
1654  s=0.d0
1655  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1656  if(lvc(idir,jdir,kdir)==1)then
1657  s(ixo^s,idir)=s(ixo^s,idir)+e(ixo^s,jdir)*b(ixo^s,kdir)
1658  else if(lvc(idir,jdir,kdir)==-1)then
1659  s(ixo^s,idir)=s(ixo^s,idir)-e(ixo^s,jdir)*b(ixo^s,kdir)
1660  end if
1661  end do; end do; end do
1662  ! equation (9)
1663  do idir = 1, ndir
1664  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))+s(ixo^s,idir)*inv_squared_c
1665  end do
1666  end subroutine mhd_to_conserved_semirelati
1667 
1668  !> Transform conservative variables into primitive ones
1669  subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
1671  integer, intent(in) :: ixI^L, ixO^L
1672  double precision, intent(inout) :: w(ixI^S, nw)
1673  double precision, intent(in) :: x(ixI^S, 1:ndim)
1674 
1675  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1676  integer :: idir
1677 
1678  if (fix_small_values) then
1679  ! fix small values preventing NaN numbers in the following converting
1680  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_origin')
1681  end if
1682 
1683  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1684 
1685  if(mhd_boris_simplification) then
1686  gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1687  ! Convert momentum to velocity
1688  do idir = 1, ndir
1689  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*gamma2
1690  end do
1691  else
1692  ! Convert momentum to velocity
1693  do idir = 1, ndir
1694  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1695  end do
1696  end if
1697 
1698  ! Calculate pressure = (gamma-1) * (e-ek-eb)
1699  if(mhd_energy) then
1700  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1701  -0.5d0*w(ixo^s,rho_)*sum(w(ixo^s,mom(:))**2,dim=ndim+1)&
1702  -mhd_mag_en(w,ixi^l,ixo^l))
1703  end if
1704 
1705  end subroutine mhd_to_primitive_origin
1706 
1707  !> Transform conservative variables into primitive ones
1708  subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
1710  integer, intent(in) :: ixI^L, ixO^L
1711  double precision, intent(inout) :: w(ixI^S, nw)
1712  double precision, intent(in) :: x(ixI^S, 1:ndim)
1713 
1714  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1715  integer :: idir
1716 
1717  if (fix_small_values) then
1718  ! fix small values preventing NaN numbers in the following converting
1719  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_hde')
1720  end if
1721 
1722  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1723 
1724  if(mhd_boris_simplification) then
1725  gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1726  ! Convert momentum to velocity
1727  do idir = 1, ndir
1728  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*gamma2
1729  end do
1730  else
1731  ! Convert momentum to velocity
1732  do idir = 1, ndir
1733  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1734  end do
1735  end if
1736 
1737  ! Calculate pressure = (gamma-1) * (e-ek)
1738  if(mhd_energy) then
1739  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-0.5d0*w(ixo^s,rho_)*sum(w(ixo^s,mom(:))**2,dim=ndim+1))
1740  end if
1741 
1742  end subroutine mhd_to_primitive_hde
1743 
1744  !> Transform conservative variables into primitive ones
1745  subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
1747  integer, intent(in) :: ixI^L, ixO^L
1748  double precision, intent(inout) :: w(ixI^S, nw)
1749  double precision, intent(in) :: x(ixI^S, 1:ndim)
1750 
1751  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1752  integer :: idir
1753 
1754  if (fix_small_values) then
1755  ! fix small values preventing NaN numbers in the following converting
1756  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_inte')
1757  end if
1758 
1759  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1760 
1761  ! Calculate pressure = (gamma-1) * e_internal
1762  if(mhd_energy) then
1763  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
1764  end if
1765 
1766  if(mhd_boris_simplification) then
1767  gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1768  ! Convert momentum to velocity
1769  do idir = 1, ndir
1770  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*gamma2
1771  end do
1772  else
1773  ! Convert momentum to velocity
1774  do idir = 1, ndir
1775  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1776  end do
1777  end if
1778 
1779  end subroutine mhd_to_primitive_inte
1780 
1781  !> Transform conservative variables into primitive ones
1782  subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
1784  integer, intent(in) :: ixI^L, ixO^L
1785  double precision, intent(inout) :: w(ixI^S, nw)
1786  double precision, intent(in) :: x(ixI^S, 1:ndim)
1787 
1788  double precision :: inv_rho(ixO^S), gamma2(ixO^S)
1789  integer :: idir
1790 
1791  if (fix_small_values) then
1792  ! fix small values preventing NaN numbers in the following converting
1793  call mhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_split_rho')
1794  end if
1795 
1796  inv_rho(ixo^s) = 1d0/(w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i))
1797 
1798  if(mhd_boris_simplification) then
1799  gamma2=inv_rho/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*inv_rho*inv_squared_c)
1800  ! Convert momentum to velocity
1801  do idir = 1, ndir
1802  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*gamma2
1803  end do
1804  else
1805  ! Convert momentum to velocity
1806  do idir = 1, ndir
1807  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))*inv_rho
1808  end do
1809  end if
1810 
1811  ! Calculate pressure = (gamma-1) * (e-ek-eb)
1812  if(mhd_energy) then
1813  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1814  -0.5d0*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))*&
1815  sum(w(ixo^s,mom(:))**2,dim=ndim+1)&
1816  -mhd_mag_en(w,ixi^l,ixo^l))
1817  end if
1818 
1819  end subroutine mhd_to_primitive_split_rho
1820 
1821  !> Transform conservative variables into primitive ones
1822  subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
1824  integer, intent(in) :: ixI^L, ixO^L
1825  double precision, intent(inout) :: w(ixI^S, nw)
1826  double precision, intent(in) :: x(ixI^S, 1:ndim)
1827 
1828  double precision :: inv_rho(ixO^S)
1829  double precision :: b(ixO^S,1:ndir), Ba(ixO^S,1:ndir),tmp(ixO^S), b2(ixO^S), gamma2(ixO^S)
1830  integer :: idir, jdir, kdir
1831 
1832  if (fix_small_values) then
1833  ! fix small values preventing NaN numbers in the following converting
1834  call mhd_handle_small_values_semirelati(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_semirelati')
1835  end if
1836 
1837  if(b0field) then
1838  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
1839  else
1840  ba(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
1841  end if
1842 
1843  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
1844 
1845  b2(ixo^s)=sum(ba(ixo^s,:)**2,dim=ndim+1)
1846  tmp(ixo^s)=sqrt(b2(ixo^s))
1847  where(tmp(ixo^s)>smalldouble)
1848  tmp(ixo^s)=1.d0/tmp(ixo^s)
1849  else where
1850  tmp(ixo^s)=0.d0
1851  end where
1852  do idir=1,ndir
1853  b(ixo^s,idir)=ba(ixo^s,idir)*tmp(ixo^s)
1854  end do
1855  tmp(ixo^s)=sum(b(ixo^s,:)*w(ixo^s,mom(:)),dim=ndim+1)
1856 
1857  ! Va^2/c^2
1858  b2(ixo^s)=b2(ixo^s)*inv_rho(ixo^s)*inv_squared_c
1859  ! equation (15)
1860  gamma2(ixo^s)=1.d0/(1.d0+b2(ixo^s))
1861  ! Convert momentum to velocity
1862  do idir = 1, ndir
1863  w(ixo^s, mom(idir)) = gamma2*(w(ixo^s, mom(idir))+b2*b(ixo^s,idir)*tmp)*inv_rho
1864  end do
1865 
1866  if(mhd_internal_e) then
1867  ! internal energy to pressure
1868  w(ixo^s,p_)=gamma_1*w(ixo^s,e_)
1869  else if(mhd_energy) then
1870  ! E=Bxv
1871  b=0.d0
1872  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
1873  if(lvc(idir,jdir,kdir)==1)then
1874  b(ixo^s,idir)=b(ixo^s,idir)+ba(ixo^s,jdir)*w(ixo^s,mom(kdir))
1875  else if(lvc(idir,jdir,kdir)==-1)then
1876  b(ixo^s,idir)=b(ixo^s,idir)-ba(ixo^s,jdir)*w(ixo^s,mom(kdir))
1877  end if
1878  end do; end do; end do
1879  ! Calculate pressure = (gamma-1) * (e-eK-eB-eE)
1880  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
1881  -half*(sum(w(ixo^s,mom(:))**2,dim=ndim+1)*w(ixo^s,rho_)&
1882  +sum(w(ixo^s,mag(:))**2,dim=ndim+1)&
1883  +sum(b(ixo^s,:)**2,dim=ndim+1)*inv_squared_c))
1884  end if
1885 
1886  end subroutine mhd_to_primitive_semirelati
1887 
1888  !> Transform internal energy to total energy
1889  subroutine mhd_ei_to_e(ixI^L,ixO^L,w,x)
1891  integer, intent(in) :: ixi^l, ixo^l
1892  double precision, intent(inout) :: w(ixi^s, nw)
1893  double precision, intent(in) :: x(ixi^s, 1:ndim)
1894 
1895  ! Calculate total energy from internal, kinetic and magnetic energy
1896  w(ixi^s,e_)=w(ixi^s,e_)&
1897  +mhd_kin_en(w,ixi^l,ixi^l)&
1898  +mhd_mag_en(w,ixi^l,ixi^l)
1899 
1900  end subroutine mhd_ei_to_e
1901 
1902  !> Transform internal energy to hydrodynamic energy
1903  subroutine mhd_ei_to_e_hde(ixI^L,ixO^L,w,x)
1905  integer, intent(in) :: ixI^L, ixO^L
1906  double precision, intent(inout) :: w(ixI^S, nw)
1907  double precision, intent(in) :: x(ixI^S, 1:ndim)
1908 
1909  ! Calculate hydrodynamic energy from internal and kinetic
1910  w(ixi^s,e_)=w(ixi^s,e_)+mhd_kin_en(w,ixi^l,ixi^l)
1911 
1912  end subroutine mhd_ei_to_e_hde
1913 
1914  !> Transform internal energy to total energy and velocity to momentum
1915  subroutine mhd_ei_to_e_semirelati(ixI^L,ixO^L,w,x)
1917  integer, intent(in) :: ixI^L, ixO^L
1918  double precision, intent(inout) :: w(ixI^S, nw)
1919  double precision, intent(in) :: x(ixI^S, 1:ndim)
1920 
1921  w(ixi^s,p_)=w(ixi^s,e_)*gamma_1
1922  call mhd_to_conserved_semirelati(ixi^l,ixi^l,w,x)
1923 
1924  end subroutine mhd_ei_to_e_semirelati
1925 
1926  !> Transform total energy to internal energy
1927  subroutine mhd_e_to_ei(ixI^L,ixO^L,w,x)
1929  integer, intent(in) :: ixi^l, ixo^l
1930  double precision, intent(inout) :: w(ixi^s, nw)
1931  double precision, intent(in) :: x(ixi^s, 1:ndim)
1932 
1933  ! Calculate ei = e - ek - eb
1934  w(ixi^s,e_)=w(ixi^s,e_)&
1935  -mhd_kin_en(w,ixi^l,ixi^l)&
1936  -mhd_mag_en(w,ixi^l,ixi^l)
1937 
1938  if(fix_small_values) then
1939  call mhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'mhd_e_to_ei')
1940  end if
1941 
1942  end subroutine mhd_e_to_ei
1943 
1944  !> Transform hydrodynamic energy to internal energy
1945  subroutine mhd_e_to_ei_hde(ixI^L,ixO^L,w,x)
1947  integer, intent(in) :: ixI^L, ixO^L
1948  double precision, intent(inout) :: w(ixI^S, nw)
1949  double precision, intent(in) :: x(ixI^S, 1:ndim)
1950 
1951  ! Calculate ei = e - ek
1952  w(ixi^s,e_)=w(ixi^s,e_)-mhd_kin_en(w,ixi^l,ixi^l)
1953 
1954  if(fix_small_values) then
1955  call mhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'mhd_e_to_ei_hde')
1956  end if
1957 
1958  end subroutine mhd_e_to_ei_hde
1959 
1960  !> Transform total energy to internal energy and momentum to velocity
1961  subroutine mhd_e_to_ei_semirelati(ixI^L,ixO^L,w,x)
1963  integer, intent(in) :: ixI^L, ixO^L
1964  double precision, intent(inout) :: w(ixI^S, nw)
1965  double precision, intent(in) :: x(ixI^S, 1:ndim)
1966 
1967  call mhd_to_primitive_semirelati(ixi^l,ixi^l,w,x)
1968  w(ixi^s,e_)=w(ixi^s,p_)*inv_gamma_1
1969 
1970  end subroutine mhd_e_to_ei_semirelati
1971 
1972  subroutine mhd_handle_small_values_semirelati(primitive, w, x, ixI^L, ixO^L, subname)
1974  use mod_small_values
1975  logical, intent(in) :: primitive
1976  integer, intent(in) :: ixI^L,ixO^L
1977  double precision, intent(inout) :: w(ixI^S,1:nw)
1978  double precision, intent(in) :: x(ixI^S,1:ndim)
1979  character(len=*), intent(in) :: subname
1980 
1981  double precision :: pressure(ixI^S), inv_rho(ixI^S), b2(ixI^S), tmp(ixI^S), gamma2(ixI^S)
1982  double precision :: b(ixI^S,1:ndir), v(ixI^S,1:ndir), Ba(ixI^S,1:ndir)
1983  integer :: idir, jdir, kdir, ix^D
1984  logical :: flag(ixI^S,1:nw)
1985 
1986  flag=.false.
1987  where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
1988 
1989  if(mhd_energy) then
1990  if(primitive) then
1991  where(w(ixo^s,p_) < small_pressure) flag(ixo^s,e_) = .true.
1992  else
1993  if(mhd_internal_e) then
1994  pressure(ixi^s)=gamma_1*w(ixi^s,e_)
1995  where(pressure(ixo^s) < small_pressure) flag(ixo^s,p_) = .true.
1996  else
1997  if(b0field) then
1998  ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,b0i)
1999  else
2000  ba(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
2001  end if
2002  inv_rho(ixi^s) = 1d0/w(ixi^s,rho_)
2003  b2(ixi^s)=sum(ba(ixi^s,:)**2,dim=ndim+1)
2004  tmp(ixi^s)=sqrt(b2(ixi^s))
2005  where(tmp(ixi^s)>smalldouble)
2006  tmp(ixi^s)=1.d0/tmp(ixi^s)
2007  else where
2008  tmp(ixi^s)=0.d0
2009  end where
2010  do idir=1,ndir
2011  b(ixi^s,idir)=ba(ixi^s,idir)*tmp(ixi^s)
2012  end do
2013  tmp(ixi^s)=sum(b(ixi^s,:)*w(ixi^s,mom(:)),dim=ndim+1)
2014  ! Va^2/c^2
2015  b2(ixi^s)=b2(ixi^s)*inv_rho(ixi^s)*inv_squared_c
2016  ! equation (15)
2017  gamma2(ixi^s)=1.d0/(1.d0+b2(ixi^s))
2018  ! Convert momentum to velocity
2019  do idir = 1, ndir
2020  v(ixi^s,idir) = gamma2*(w(ixi^s, mom(idir))+b2*b(ixi^s,idir)*tmp(ixi^s))*inv_rho(ixi^s)
2021  end do
2022  ! E=Bxv
2023  b=0.d0
2024  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
2025  if(lvc(idir,jdir,kdir)==1)then
2026  b(ixi^s,idir)=b(ixi^s,idir)+ba(ixi^s,jdir)*v(ixi^s,kdir)
2027  else if(lvc(idir,jdir,kdir)==-1)then
2028  b(ixi^s,idir)=b(ixi^s,idir)-ba(ixi^s,jdir)*v(ixi^s,kdir)
2029  end if
2030  end do; end do; end do
2031  ! Calculate pressure p = (gamma-1)(e-eK-eB-eE)
2032  pressure(ixi^s)=gamma_1*(w(ixi^s,e_)&
2033  -half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,rho_)&
2034  +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
2035  +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c))
2036  where(pressure(ixo^s) < small_pressure) flag(ixo^s,p_) = .true.
2037  end if
2038  end if
2039  end if
2040 
2041  if(any(flag)) then
2042  select case (small_values_method)
2043  case ("replace")
2044  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2045  if(mhd_energy) then
2046  if(primitive) then
2047  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
2048  else
2049  if(mhd_internal_e) then
2050  ! internal energy
2051  {do ix^db=ixomin^db,ixomax^db\}
2052  if(flag(ix^d,e_)) then
2053  w(ix^d,e_)=small_pressure*inv_gamma_1
2054  end if
2055  {end do\}
2056  else
2057  {do ix^db=ixomin^db,ixomax^db\}
2058  if(flag(ix^d,e_)) then
2059  w(ix^d,e_)=small_pressure*inv_gamma_1+half*(sum(v(ix^d,:)**2)*w(ix^d,rho_)&
2060  +sum(w(ix^d,mag(:))**2)+sum(b(ix^d,:)**2)*inv_squared_c)
2061  end if
2062  {end do\}
2063  end if
2064  end if
2065  end if
2066  case ("average")
2067  ! do averaging of density
2068  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2069  if(mhd_energy) then
2070  if(primitive) then
2071  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2072  else
2073  if(mhd_internal_e) then
2074  ! internal energy
2075  w(ixi^s,e_)=pressure(ixi^s)
2076  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2077  w(ixi^s,e_)=w(ixi^s,p_)*inv_gamma_1
2078  else
2079  w(ixi^s,e_)=pressure(ixi^s)
2080  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2081  w(ixi^s,e_)=w(ixi^s,p_)*inv_gamma_1&
2082  +half*(sum(v(ixi^s,:)**2,dim=ndim+1)*w(ixi^s,rho_)&
2083  +sum(w(ixi^s,mag(:))**2,dim=ndim+1)&
2084  +sum(b(ixi^s,:)**2,dim=ndim+1)*inv_squared_c)
2085  end if
2086  end if
2087  end if
2088  case default
2089  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2090  end select
2091  end if
2092  end subroutine mhd_handle_small_values_semirelati
2093 
2094  subroutine mhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
2096  use mod_small_values
2097  logical, intent(in) :: primitive
2098  integer, intent(in) :: ixI^L,ixO^L
2099  double precision, intent(inout) :: w(ixI^S,1:nw)
2100  double precision, intent(in) :: x(ixI^S,1:ndim)
2101  character(len=*), intent(in) :: subname
2102 
2103  integer :: idir
2104  logical :: flag(ixI^S,1:nw)
2105  double precision :: tmp2(ixI^S)
2106 
2107  call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2108 
2109  if(any(flag)) then
2110  select case (small_values_method)
2111  case ("replace")
2112  if(has_equi_rho0) then
2113  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = &
2114  small_density-block%equi_vars(ixo^s,equi_rho0_,0)
2115  else
2116  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2117  end if
2118  do idir = 1, ndir
2119  if(small_values_fix_iw(mom(idir))) then
2120  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
2121  end if
2122  end do
2123  if(mhd_energy) then
2124  if(primitive) then
2125  if(has_equi_pe0) then
2126  tmp2(ixo^s) = small_pressure - &
2127  block%equi_vars(ixo^s,equi_pe0_,0)
2128  where(flag(ixo^s,e_)) w(ixo^s,p_) = tmp2(ixo^s)
2129  else
2130  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
2131  end if
2132  else
2133  ! conserved
2134  if(has_equi_pe0) then
2135  tmp2(ixo^s) = small_e - &
2136  block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
2137  where(flag(ixo^s,e_))
2138  w(ixo^s,e_) = tmp2(ixo^s)+&
2139  mhd_kin_en(w,ixi^l,ixo^l)+&
2140  mhd_mag_en(w,ixi^l,ixo^l)
2141  end where
2142  else
2143  where(flag(ixo^s,e_))
2144  w(ixo^s,e_) = small_e+&
2145  mhd_kin_en(w,ixi^l,ixo^l)+&
2146  mhd_mag_en(w,ixi^l,ixo^l)
2147  end where
2148  end if
2149  end if
2150  end if
2151  case ("average")
2152  ! do averaging of density
2153  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2154  if(mhd_energy) then
2155  if(primitive)then
2156  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2157  else
2158  ! do averaging of internal energy
2159  w(ixi^s,e_)=w(ixi^s,e_)&
2160  -mhd_kin_en(w,ixi^l,ixi^l)&
2161  -mhd_mag_en(w,ixi^l,ixi^l)
2162  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
2163  ! convert back
2164  w(ixi^s,e_)=w(ixi^s,e_)&
2165  +mhd_kin_en(w,ixi^l,ixi^l)&
2166  +mhd_mag_en(w,ixi^l,ixi^l)
2167  end if
2168  end if
2169  case default
2170  if(.not.primitive) then
2171  !convert w to primitive
2172  ! Calculate pressure = (gamma-1) * (e-ek-eb)
2173  if(mhd_energy) then
2174  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)&
2175  -mhd_kin_en(w,ixi^l,ixo^l)&
2176  -mhd_mag_en(w,ixi^l,ixo^l))
2177  end if
2178  ! Convert momentum to velocity
2179  if(has_equi_rho0) then
2180  tmp2(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
2181  else
2182  tmp2(ixo^s) = w(ixo^s,rho_)
2183  end if
2184  do idir = 1, ndir
2185  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/tmp2(ixo^s)
2186  end do
2187  end if
2188  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2189  end select
2190  end if
2191 
2192  end subroutine mhd_handle_small_values_origin
2193 
2194  subroutine mhd_handle_small_values_inte(primitive, w, x, ixI^L, ixO^L, subname)
2196  use mod_small_values
2197  logical, intent(in) :: primitive
2198  integer, intent(in) :: ixI^L,ixO^L
2199  double precision, intent(inout) :: w(ixI^S,1:nw)
2200  double precision, intent(in) :: x(ixI^S,1:ndim)
2201  character(len=*), intent(in) :: subname
2202 
2203  integer :: idir
2204  logical :: flag(ixI^S,1:nw)
2205  double precision :: tmp2(ixI^S)
2206 
2207  call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2208 
2209  if(any(flag)) then
2210  select case (small_values_method)
2211  case ("replace")
2212  if(has_equi_rho0) then
2213  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = &
2214  small_density-block%equi_vars(ixo^s,equi_rho0_,0)
2215  else
2216  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2217  end if
2218  do idir = 1, ndir
2219  if(small_values_fix_iw(mom(idir))) then
2220  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
2221  end if
2222  end do
2223  if(mhd_energy) then
2224  if(primitive) then
2225  if(has_equi_pe0) then
2226  tmp2(ixo^s) = small_pressure - &
2227  block%equi_vars(ixo^s,equi_pe0_,0)
2228  where(flag(ixo^s,e_)) w(ixo^s,p_) = tmp2(ixo^s)
2229  else
2230  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
2231  end if
2232  else
2233  ! conserved
2234  if(has_equi_pe0) then
2235  tmp2(ixo^s) = small_e - &
2236  block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
2237  where(flag(ixo^s,e_))
2238  w(ixo^s,e_)=tmp2(ixo^s)
2239  end where
2240  else
2241  where(flag(ixo^s,e_))
2242  w(ixo^s,e_)=small_e
2243  end where
2244  end if
2245  end if
2246  end if
2247  case ("average")
2248  ! do averaging of density
2249  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2250  if(mhd_energy) then
2251  if(primitive)then
2252  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2253  else
2254  ! do averaging of internal energy
2255  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
2256  end if
2257  end if
2258  case default
2259  if(.not.primitive) then
2260  !convert w to primitive
2261  if(mhd_energy) then
2262  w(ixo^s,p_)=w(ixo^s,e_)*gamma_1
2263  end if
2264  ! Convert momentum to velocity
2265  if(has_equi_rho0) then
2266  tmp2(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,0)
2267  else
2268  tmp2(ixo^s) = w(ixo^s,rho_)
2269  end if
2270  do idir = 1, ndir
2271  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/tmp2(ixo^s)
2272  end do
2273  end if
2274  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2275  end select
2276  end if
2277 
2278  end subroutine mhd_handle_small_values_inte
2279 
2280  subroutine mhd_handle_small_values_hde(primitive, w, x, ixI^L, ixO^L, subname)
2282  use mod_small_values
2283  logical, intent(in) :: primitive
2284  integer, intent(in) :: ixI^L,ixO^L
2285  double precision, intent(inout) :: w(ixI^S,1:nw)
2286  double precision, intent(in) :: x(ixI^S,1:ndim)
2287  character(len=*), intent(in) :: subname
2288 
2289  integer :: idir
2290  logical :: flag(ixI^S,1:nw)
2291  double precision :: tmp2(ixI^S)
2292 
2293  call phys_check_w(primitive, ixi^l, ixo^l, w, flag)
2294 
2295  if(any(flag)) then
2296  select case (small_values_method)
2297  case ("replace")
2298  where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
2299  do idir = 1, ndir
2300  if(small_values_fix_iw(mom(idir))) then
2301  where(flag(ixo^s,rho_)) w(ixo^s, mom(idir)) = 0.0d0
2302  end if
2303  end do
2304 
2305  if(mhd_energy) then
2306  if(primitive) then
2307  where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
2308  else
2309  where(flag(ixo^s,e_))
2310  w(ixo^s,e_) = small_e+mhd_kin_en(w,ixi^l,ixo^l)
2311  end where
2312  end if
2313  end if
2314  case ("average")
2315  ! do averaging of density
2316  call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
2317  if(mhd_energy) then
2318  if(primitive) then
2319  call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
2320  else
2321  ! do averaging of internal energy
2322  w(ixi^s,e_)=w(ixi^s,e_)-mhd_kin_en(w,ixi^l,ixi^l)
2323  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
2324  ! convert back
2325  w(ixi^s,e_)=w(ixi^s,e_)+mhd_kin_en(w,ixi^l,ixi^l)
2326  end if
2327  end if
2328  case default
2329  if(.not.primitive) then
2330  !convert w to primitive
2331  ! Calculate pressure = (gamma-1) * (e-ek)
2332  if(mhd_energy) then
2333  w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l))
2334  end if
2335  ! Convert momentum to velocity
2336  do idir = 1, ndir
2337  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/w(ixo^s,rho_)
2338  end do
2339  end if
2340  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
2341  end select
2342  end if
2343 
2344  end subroutine mhd_handle_small_values_hde
2345 
2346  !> Calculate v vector
2347  subroutine mhd_get_v_origin(w,x,ixI^L,ixO^L,v)
2349 
2350  integer, intent(in) :: ixI^L, ixO^L
2351  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2352  double precision, intent(out) :: v(ixI^S,ndir)
2353 
2354  double precision :: rho(ixI^S)
2355  integer :: idir
2356 
2357  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2358 
2359  rho(ixo^s)=1.d0/rho(ixo^s)
2360  ! Convert momentum to velocity
2361  do idir = 1, ndir
2362  v(ixo^s, idir) = w(ixo^s, mom(idir))*rho(ixo^s)
2363  end do
2364 
2365  end subroutine mhd_get_v_origin
2366 
2367  !> Calculate v vector
2368  subroutine mhd_get_v_boris(w,x,ixI^L,ixO^L,v)
2370 
2371  integer, intent(in) :: ixI^L, ixO^L
2372  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
2373  double precision, intent(out) :: v(ixI^S,ndir)
2374 
2375  double precision :: rho(ixI^S), gamma2(ixO^S)
2376  integer :: idir
2377 
2378  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2379 
2380  rho(ixo^s)=1.d0/rho(ixo^s)
2381  gamma2=1.d0/(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)*rho(ixo^s)*inv_squared_c)
2382  ! Convert momentum to velocity
2383  do idir = 1, ndir
2384  v(ixo^s, idir) = w(ixo^s, mom(idir))*rho(ixo^s)*gamma2
2385  end do
2386 
2387  end subroutine mhd_get_v_boris
2388 
2389  !> Calculate v component
2390  subroutine mhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
2392 
2393  integer, intent(in) :: ixi^l, ixo^l, idim
2394  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
2395  double precision, intent(out) :: v(ixi^s)
2396 
2397  double precision :: rho(ixi^s)
2398 
2399  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
2400 
2401  if(mhd_boris_simplification) then
2402  v(ixo^s) = w(ixo^s, mom(idim)) / rho(ixo^s) &
2403  /(1.d0+sum(w(ixo^s,mag(:))**2,dim=ndim+1)/rho(ixo^s)*inv_squared_c)
2404  else
2405  v(ixo^s) = w(ixo^s, mom(idim)) / rho(ixo^s)
2406  end if
2407 
2408  end subroutine mhd_get_v_idim
2409 
2410  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
2411  subroutine mhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
2413 
2414  integer, intent(in) :: ixI^L, ixO^L, idim
2415  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2416  double precision, intent(inout) :: cmax(ixI^S)
2417  double precision :: vel(ixI^S)
2418 
2419  call mhd_get_csound(w,x,ixi^l,ixo^l,idim,cmax)
2420  call mhd_get_v_idim(w,x,ixi^l,ixo^l,idim,vel)
2421 
2422  cmax(ixo^s)=abs(vel(ixo^s))+cmax(ixo^s)
2423 
2424  end subroutine mhd_get_cmax_origin
2425 
2426  !> Calculate cmax_idim for semirelativistic MHD
2427  subroutine mhd_get_cmax_semirelati(w,x,ixI^L,ixO^L,idim,cmax)
2429 
2430  integer, intent(in) :: ixI^L, ixO^L, idim
2431  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2432  double precision, intent(inout):: cmax(ixI^S)
2433  double precision :: wprim(ixI^S,nw)
2434  double precision :: csound(ixO^S), AvMinCs2(ixO^S), idim_Alfven_speed2(ixO^S)
2435  double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), gamma2(ixO^S), B(ixO^S,1:ndir)
2436 
2437  if(b0field) then
2438  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
2439  else
2440  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
2441  end if
2442  inv_rho = 1.d0/w(ixo^s,rho_)
2443 
2444  alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
2445  gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
2446 
2447  wprim=w
2448  call mhd_to_primitive(ixi^l,ixo^l,wprim,x)
2449  cmax(ixo^s)=1.d0-gamma2*wprim(ixo^s,mom(idim))**2*inv_squared_c
2450  ! equation (69)
2451  alfven_speed2=alfven_speed2*cmax(ixo^s)
2452 
2453  ! squared sound speed
2454  if(mhd_energy) then
2455  csound=mhd_gamma*wprim(ixo^s,p_)*inv_rho
2456  else
2457  csound=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
2458  end if
2459 
2460  idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
2461 
2462  ! Va_hat^2+a_hat^2 equation (57)
2463  alfven_speed2=alfven_speed2+csound*(1.d0+idim_alfven_speed2*inv_squared_c)
2464 
2465  avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound*idim_alfven_speed2*cmax(ixo^s)
2466 
2467  where(avmincs2<zero)
2468  avmincs2=zero
2469  end where
2470 
2471  ! equation (68) fast magnetosonic wave speed
2472  csound = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
2473  cmax(ixo^s)=gamma2*abs(wprim(ixo^s,mom(idim)))+csound
2474 
2475  end subroutine mhd_get_cmax_semirelati
2476 
2477  subroutine mhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
2479 
2480  integer, intent(in) :: ixI^L, ixO^L
2481  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2482  double precision, intent(inout) :: a2max(ndim)
2483  double precision :: a2(ixI^S,ndim,nw)
2484  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
2485 
2486  a2=zero
2487  do i = 1,ndim
2488  !> 4th order
2489  hxo^l=ixo^l-kr(i,^d);
2490  gxo^l=hxo^l-kr(i,^d);
2491  jxo^l=ixo^l+kr(i,^d);
2492  kxo^l=jxo^l+kr(i,^d);
2493  a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
2494  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
2495  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
2496  end do
2497  end subroutine mhd_get_a2max
2498 
2499  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
2500  subroutine mhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
2502  use mod_geometry
2503  integer, intent(in) :: ixI^L,ixO^L
2504  double precision, intent(in) :: x(ixI^S,1:ndim)
2505  double precision, intent(inout) :: w(ixI^S,1:nw)
2506  double precision, intent(out) :: Tco_local,Tmax_local
2507 
2508  double precision, parameter :: trac_delta=0.25d0
2509  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
2510  double precision, dimension(ixI^S,1:ndir) :: bunitvec
2511  double precision, dimension(ixI^S,1:ndim) :: gradT
2512  double precision :: Bdir(ndim)
2513  double precision :: ltrc,ltrp,altr(ixI^S)
2514  integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
2515  integer :: jxP^L,hxP^L,ixP^L,ixQ^L
2516  logical :: lrlt(ixI^S)
2517 
2518  call mhd_get_temperature(w,x,ixi^l,ixi^l,te)
2519  tco_local=zero
2520  tmax_local=maxval(te(ixo^s))
2521 
2522  {^ifoned
2523  select case(mhd_trac_type)
2524  case(0)
2525  !> test case, fixed cutoff temperature
2526  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
2527  case(1)
2528  hxo^l=ixo^l-1;
2529  jxo^l=ixo^l+1;
2530  lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
2531  lrlt=.false.
2532  where(lts(ixo^s) > trac_delta)
2533  lrlt(ixo^s)=.true.
2534  end where
2535  if(any(lrlt(ixo^s))) then
2536  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
2537  end if
2538  case(2)
2539  !> iijima et al. 2021, LTRAC method
2540  ltrc=1.5d0
2541  ltrp=4.d0
2542  ixp^l=ixo^l^ladd1;
2543  hxo^l=ixo^l-1;
2544  jxo^l=ixo^l+1;
2545  hxp^l=ixp^l-1;
2546  jxp^l=ixp^l+1;
2547  lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
2548  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2549  lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
2550  block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
2551  case default
2552  call mpistop("mhd_trac_type not allowed for 1D simulation")
2553  end select
2554  }
2555  {^nooned
2556  select case(mhd_trac_type)
2557  case(0)
2558  !> test case, fixed cutoff temperature
2559  block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
2560  case(1,4,6)
2561  ! temperature gradient at cell centers
2562  do idims=1,ndim
2563  call gradient(te,ixi^l,ixo^l,idims,tmp1)
2564  gradt(ixo^s,idims)=tmp1(ixo^s)
2565  end do
2566  ! B vector
2567  if(b0field) then
2568  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
2569  else
2570  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
2571  end if
2572  if(mhd_trac_type .gt. 1) then
2573  ! B direction at cell center
2574  bdir=zero
2575  {do ixa^d=0,1\}
2576  ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
2577  bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
2578  {end do\}
2579  if(sum(bdir(:)**2) .gt. zero) then
2580  bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
2581  end if
2582  block%special_values(3:ndim+2)=bdir(1:ndim)
2583  end if
2584  tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
2585  where(tmp1(ixo^s)/=0.d0)
2586  tmp1(ixo^s)=1.d0/tmp1(ixo^s)
2587  elsewhere
2588  tmp1(ixo^s)=bigdouble
2589  end where
2590  ! b unit vector: magnetic field direction vector
2591  do idims=1,ndim
2592  bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
2593  end do
2594  ! temperature length scale inversed
2595  lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
2596  ! fraction of cells size to temperature length scale
2597  if(slab_uniform) then
2598  lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
2599  else
2600  lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
2601  end if
2602  lrlt=.false.
2603  where(lts(ixo^s) > trac_delta)
2604  lrlt(ixo^s)=.true.
2605  end where
2606  if(any(lrlt(ixo^s))) then
2607  block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
2608  else
2609  block%special_values(1)=zero
2610  end if
2611  block%special_values(2)=tmax_local
2612  case(2)
2613  !> iijima et al. 2021, LTRAC method
2614  ltrc=1.5d0
2615  ltrp=4.d0
2616  ixp^l=ixo^l^ladd2;
2617  ! temperature gradient at cell centers
2618  do idims=1,ndim
2619  ixq^l=ixp^l;
2620  hxp^l=ixp^l;
2621  jxp^l=ixp^l;
2622  select case(idims)
2623  {case(^d)
2624  ixqmin^d=ixqmin^d+1
2625  ixqmax^d=ixqmax^d-1
2626  hxpmax^d=ixpmin^d
2627  jxpmin^d=ixpmax^d
2628  \}
2629  end select
2630  call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
2631  call gradientx(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),.false.)
2632  call gradientq(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims))
2633  end do
2634  ! B vector
2635  if(b0field) then
2636  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
2637  else
2638  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
2639  end if
2640  tmp1(ixp^s)=1.d0/(dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))+smalldouble)
2641  ! b unit vector: magnetic field direction vector
2642  do idims=1,ndim
2643  bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
2644  end do
2645  ! temperature length scale inversed
2646  lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
2647  ! fraction of cells size to temperature length scale
2648  if(slab_uniform) then
2649  lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
2650  else
2651  lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
2652  end if
2653  lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
2654 
2655  altr=zero
2656  ! need one ghost layer for thermal conductivity
2657  ixp^l=ixo^l^ladd1;
2658  do idims=1,ndim
2659  hxo^l=ixp^l-kr(idims,^d);
2660  jxo^l=ixp^l+kr(idims,^d);
2661  altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
2662  end do
2663  block%wextra(ixp^s,tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
2664  case(3,5)
2665  !> do nothing here
2666  case default
2667  call mpistop("unknown mhd_trac_type")
2668  end select
2669  }
2670  end subroutine mhd_get_tcutoff
2671 
2672  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
2673  subroutine mhd_get_h_speed(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2675 
2676  integer, intent(in) :: ixI^L, ixO^L, idim
2677  double precision, intent(in) :: wprim(ixI^S, nw)
2678  double precision, intent(in) :: x(ixI^S,1:ndim)
2679  double precision, intent(out) :: Hspeed(ixI^S,1:number_species)
2680 
2681  double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2682  integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2683 
2684  hspeed=0.d0
2685  ixa^l=ixo^l^ladd1;
2686  do id=1,ndim
2687  call mhd_get_csound_prim(wprim,x,ixi^l,ixa^l,id,tmp)
2688  csound(ixa^s,id)=tmp(ixa^s)
2689  end do
2690  ixcmax^d=ixomax^d;
2691  ixcmin^d=ixomin^d+kr(idim,^d)-1;
2692  jxcmax^d=ixcmax^d+kr(idim,^d);
2693  jxcmin^d=ixcmin^d+kr(idim,^d);
2694  hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,mom(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom(idim))+csound(ixc^s,idim))
2695 
2696  do id=1,ndim
2697  if(id==idim) cycle
2698  ixamax^d=ixcmax^d+kr(id,^d);
2699  ixamin^d=ixcmin^d+kr(id,^d);
2700  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(ixc^s,mom(id))+csound(ixc^s,id)))
2701  ixamax^d=ixcmax^d-kr(id,^d);
2702  ixamin^d=ixcmin^d-kr(id,^d);
2703  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,mom(id))+csound(ixc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
2704  end do
2705 
2706  do id=1,ndim
2707  if(id==idim) cycle
2708  ixamax^d=jxcmax^d+kr(id,^d);
2709  ixamin^d=jxcmin^d+kr(id,^d);
2710  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom(id))+csound(ixa^s,id)-wprim(jxc^s,mom(id))+csound(jxc^s,id)))
2711  ixamax^d=jxcmax^d-kr(id,^d);
2712  ixamin^d=jxcmin^d-kr(id,^d);
2713  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,mom(id))+csound(jxc^s,id)-wprim(ixa^s,mom(id))+csound(ixa^s,id)))
2714  end do
2715 
2716  end subroutine mhd_get_h_speed
2717 
2718  !> Estimating bounds for the minimum and maximum signal velocities without split
2719  subroutine mhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2721 
2722  integer, intent(in) :: ixI^L, ixO^L, idim
2723  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2724  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2725  double precision, intent(in) :: x(ixI^S,1:ndim)
2726  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2727  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2728  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2729 
2730  double precision :: wmean(ixI^S,nw)
2731  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2732  integer :: ix^D
2733 
2734  select case (boundspeed)
2735  case (1)
2736  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2737  ! Methods for Fluid Dynamics" by Toro.
2738  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
2739  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
2740  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2741  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2742  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2743  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2744  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2745  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2746  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
2747  dmean(ixo^s)=sqrt(dmean(ixo^s))
2748  if(present(cmin)) then
2749  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2750  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2751  if(h_correction) then
2752  {do ix^db=ixomin^db,ixomax^db\}
2753  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2754  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2755  {end do\}
2756  end if
2757  else
2758  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2759  end if
2760  case (2)
2761  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2762  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/wmean(ixo^s,rho_)
2763  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2764  if(present(cmin)) then
2765  cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2766  cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2767  if(h_correction) then
2768  {do ix^db=ixomin^db,ixomax^db\}
2769  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2770  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2771  {end do\}
2772  end if
2773  else
2774  cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2775  end if
2776  case (3)
2777  ! Miyoshi 2005 JCP 208, 315 equation (67)
2778  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2779  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2780  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2781  if(present(cmin)) then
2782  cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2783  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2784  if(h_correction) then
2785  {do ix^db=ixomin^db,ixomax^db\}
2786  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2787  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2788  {end do\}
2789  end if
2790  else
2791  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2792  end if
2793  end select
2794 
2795  end subroutine mhd_get_cbounds
2796 
2797  !> Estimating bounds for the minimum and maximum signal velocities without split
2798  subroutine mhd_get_cbounds_semirelati(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2800 
2801  integer, intent(in) :: ixI^L, ixO^L, idim
2802  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2803  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2804  double precision, intent(in) :: x(ixI^S,1:ndim)
2805  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2806  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2807  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2808 
2809  double precision, dimension(ixO^S) :: csoundL, csoundR, gamma2L, gamma2R
2810 
2811  ! Miyoshi 2005 JCP 208, 315 equation (67)
2812  call mhd_get_csound_semirelati(wlp,x,ixi^l,ixo^l,idim,csoundl,gamma2l)
2813  call mhd_get_csound_semirelati(wrp,x,ixi^l,ixo^l,idim,csoundr,gamma2r)
2814  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2815  if(present(cmin)) then
2816  cmin(ixo^s,1)=min(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2817  cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2818  else
2819  cmax(ixo^s,1)=max(gamma2l*wlp(ixo^s,mom(idim)),gamma2r*wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2820  end if
2821 
2822  end subroutine mhd_get_cbounds_semirelati
2823 
2824  !> Estimating bounds for the minimum and maximum signal velocities with rho split
2825  subroutine mhd_get_cbounds_split_rho(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2827 
2828  integer, intent(in) :: ixI^L, ixO^L, idim
2829  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2830  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2831  double precision, intent(in) :: x(ixI^S,1:ndim)
2832  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2833  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2834  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2835 
2836  double precision :: wmean(ixI^S,nw)
2837  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2838  integer :: ix^D
2839  double precision :: rho(ixI^S)
2840 
2841  select case (boundspeed)
2842  case (1)
2843  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2844  ! Methods for Fluid Dynamics" by Toro.
2845  tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2846  tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2847  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2848  umean(ixo^s)=(wlp(ixo^s,mom(idim))*tmp1(ixo^s)+wrp(ixo^s,mom(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2849  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2850  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2851  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2852  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2853  (wrp(ixo^s,mom(idim))-wlp(ixo^s,mom(idim)))**2
2854  dmean(ixo^s)=sqrt(dmean(ixo^s))
2855  if(present(cmin)) then
2856  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2857  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2858  if(h_correction) then
2859  {do ix^db=ixomin^db,ixomax^db\}
2860  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2861  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2862  {end do\}
2863  end if
2864  else
2865  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2866  end if
2867  case (2)
2868  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2869  tmp1(ixo^s)=wmean(ixo^s,mom(idim))/(wmean(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i))
2870  call mhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2871  if(present(cmin)) then
2872  cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
2873  cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
2874  if(h_correction) then
2875  {do ix^db=ixomin^db,ixomax^db\}
2876  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2877  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2878  {end do\}
2879  end if
2880  else
2881  cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2882  end if
2883  case (3)
2884  ! Miyoshi 2005 JCP 208, 315 equation (67)
2885  call mhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2886  call mhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2887  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2888  if(present(cmin)) then
2889  cmin(ixo^s,1)=min(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))-csoundl(ixo^s)
2890  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2891  if(h_correction) then
2892  {do ix^db=ixomin^db,ixomax^db\}
2893  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2894  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2895  {end do\}
2896  end if
2897  else
2898  cmax(ixo^s,1)=max(wlp(ixo^s,mom(idim)),wrp(ixo^s,mom(idim)))+csoundl(ixo^s)
2899  end if
2900  end select
2901 
2902  end subroutine mhd_get_cbounds_split_rho
2903 
2904  !> prepare velocities for ct methods
2905  subroutine mhd_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2907 
2908  integer, intent(in) :: ixI^L, ixO^L, idim
2909  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2910  double precision, intent(in) :: cmax(ixI^S)
2911  double precision, intent(in), optional :: cmin(ixI^S)
2912  type(ct_velocity), intent(inout):: vcts
2913 
2914  integer :: idimE,idimN
2915 
2916  ! calculate velocities related to different UCT schemes
2917  select case(type_ct)
2918  case('average')
2919  case('uct_contact')
2920  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
2921  ! get average normal velocity at cell faces
2922  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
2923  case('uct_hll')
2924  if(.not.allocated(vcts%vbarC)) then
2925  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
2926  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
2927  end if
2928  ! Store magnitude of characteristics
2929  if(present(cmin)) then
2930  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2931  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2932  else
2933  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2934  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2935  end if
2936 
2937  idimn=mod(idim,ndir)+1 ! 'Next' direction
2938  idime=mod(idim+1,ndir)+1 ! Electric field direction
2939  ! Store velocities
2940  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
2941  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
2942  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2943  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2944  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2945 
2946  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
2947  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
2948  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2949  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2950  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2951  case default
2952  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
2953  end select
2954 
2955  end subroutine mhd_get_ct_velocity
2956 
2957  !> Calculate fast magnetosonic wave speed
2958  subroutine mhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2960 
2961  integer, intent(in) :: ixI^L, ixO^L, idim
2962  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2963  double precision, intent(out):: csound(ixI^S)
2964  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2965  double precision :: inv_rho(ixO^S)
2966 
2967  if(has_equi_rho0) then
2968  inv_rho(ixo^s) = 1d0/(w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i))
2969  else
2970  inv_rho(ixo^s) = 1d0/w(ixo^s,rho_)
2971  endif
2972 
2973  call mhd_get_csound2(w,x,ixi^l,ixo^l,csound)
2974 
2975  ! store |B|^2 in v
2976  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
2977 
2978  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
2979  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2980  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 * inv_rho
2981 
2982  where(avmincs2(ixo^s)<zero)
2983  avmincs2(ixo^s)=zero
2984  end where
2985 
2986  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2987 
2988  if (.not. mhd_hall) then
2989  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2990  if (mhd_boris_simplification) then
2991  ! equation (88)
2992  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
2993  end if
2994  else
2995  ! take the Hall velocity into account:
2996  ! most simple estimate, high k limit:
2997  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2998  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2999  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
3000  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
3001  end if
3002 
3003  end subroutine mhd_get_csound
3004 
3005  !> Calculate fast magnetosonic wave speed
3006  subroutine mhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
3008 
3009  integer, intent(in) :: ixI^L, ixO^L, idim
3010  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3011  double precision, intent(out):: csound(ixI^S)
3012  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
3013  double precision :: inv_rho(ixO^S)
3014  double precision :: tmp(ixI^S)
3015 
3016  call mhd_get_rho(w,x,ixi^l,ixo^l,tmp)
3017  inv_rho(ixo^s) = 1d0/tmp(ixo^s)
3018 
3019 
3020  if(mhd_energy) then
3021  if(has_equi_pe0) then
3022  csound(ixo^s) = w(ixo^s,e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)
3023  else
3024  csound(ixo^s) = w(ixo^s,e_)
3025  endif
3026  csound(ixo^s)=mhd_gamma*csound(ixo^s)*inv_rho
3027  else
3028  csound(ixo^s)=mhd_gamma*mhd_adiab*tmp(ixo^s)**gamma_1
3029  end if
3030 
3031  ! store |B|^2 in v
3032  b2(ixo^s) = mhd_mag_en_all(w,ixi^l,ixo^l)
3033  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
3034  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
3035  * mhd_mag_i_all(w,ixi^l,ixo^l,idim)**2 * inv_rho
3036 
3037  where(avmincs2(ixo^s)<zero)
3038  avmincs2(ixo^s)=zero
3039  end where
3040 
3041  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
3042 
3043  if (.not. mhd_hall) then
3044  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
3045  if (mhd_boris_simplification) then
3046  ! equation (88)
3047  csound(ixo^s) = mhd_gamma_alfven(w, ixi^l,ixo^l) * csound(ixo^s)
3048  end if
3049  else
3050  ! take the Hall velocity into account:
3051  ! most simple estimate, high k limit:
3052  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
3053  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
3054  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
3055  mhd_etah * sqrt(b2(ixo^s))*inv_rho*kmax)
3056  end if
3057 
3058  end subroutine mhd_get_csound_prim
3059 
3060  !> Calculate cmax_idim for semirelativistic MHD
3061  subroutine mhd_get_csound_semirelati(w,x,ixI^L,ixO^L,idim,csound,gamma2)
3063 
3064  integer, intent(in) :: ixI^L, ixO^L, idim
3065  ! here w is primitive variables
3066  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
3067  double precision, intent(out):: csound(ixO^S), gamma2(ixO^S)
3068  double precision :: AvMinCs2(ixO^S), kmax
3069  double precision :: inv_rho(ixO^S), Alfven_speed2(ixO^S), idim_Alfven_speed2(ixO^S),B(ixO^S,1:ndir)
3070 
3071  integer :: ix^D
3072 
3073  if(b0field) then
3074  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,b0i)
3075  else
3076  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3077  end if
3078 
3079  inv_rho = 1.d0/w(ixo^s,rho_)
3080 
3081  alfven_speed2=sum(b(ixo^s,:)**2,dim=ndim+1)*inv_rho
3082  gamma2 = 1.0d0/(1.d0+alfven_speed2*inv_squared_c)
3083 
3084  avmincs2=1.d0-gamma2*w(ixo^s,mom(idim))**2*inv_squared_c
3085  ! equatoin (69)
3086  alfven_speed2=alfven_speed2*avmincs2
3087 
3088  ! squared sound speed
3089  if(mhd_energy) then
3090  csound(ixo^s)=mhd_gamma*w(ixo^s,p_)*inv_rho
3091  else
3092  csound(ixo^s)=mhd_gamma*mhd_adiab*w(ixo^s,rho_)**gamma_1
3093  end if
3094 
3095  idim_alfven_speed2=b(ixo^s,idim)**2*inv_rho
3096 
3097  ! Va_hat^2+a_hat^2 equation (57)
3098  alfven_speed2=alfven_speed2+csound(ixo^s)*(1.d0+idim_alfven_speed2*inv_squared_c)
3099 
3100  avmincs2=(gamma2*alfven_speed2)**2-4.0d0*gamma2*csound(ixo^s)*idim_alfven_speed2*avmincs2
3101 
3102  where(avmincs2<zero)
3103  avmincs2=zero
3104  end where
3105 
3106  ! equation (68) fast magnetosonic speed
3107  csound(ixo^s) = sqrt(half*(gamma2*alfven_speed2+sqrt(avmincs2)))
3108 
3109  end subroutine mhd_get_csound_semirelati
3110 
3111  !> Calculate isothermal thermal pressure
3112  subroutine mhd_get_pthermal_iso(w,x,ixI^L,ixO^L,pth)
3114 
3115  integer, intent(in) :: ixI^L, ixO^L
3116  double precision, intent(in) :: w(ixI^S,nw)
3117  double precision, intent(in) :: x(ixI^S,1:ndim)
3118  double precision, intent(out):: pth(ixI^S)
3119 
3120  call mhd_get_rho(w,x,ixi^l,ixo^l,pth)
3121  pth(ixo^s)=mhd_adiab*pth(ixo^s)**mhd_gamma
3122 
3123  end subroutine mhd_get_pthermal_iso
3124 
3125  !> Calculate thermal pressure from internal energy
3126  subroutine mhd_get_pthermal_eint(w,x,ixI^L,ixO^L,pth)
3129 
3130  integer, intent(in) :: ixI^L, ixO^L
3131  double precision, intent(in) :: w(ixI^S,nw)
3132  double precision, intent(in) :: x(ixI^S,1:ndim)
3133  double precision, intent(out):: pth(ixI^S)
3134  integer :: iw, ix^D
3135 
3136  pth(ixo^s)=gamma_1*w(ixo^s,e_)
3137 
3138  if(has_equi_pe0) then
3139  pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe0_,b0i)
3140  end if
3141 
3142  if (fix_small_values) then
3143  {do ix^db= ixo^lim^db\}
3144  if(pth(ix^d)<small_pressure) then
3145  pth(ix^d)=small_pressure
3146  end if
3147  {enddo^d&\}
3148  else if (check_small_values) then
3149  {do ix^db= ixo^lim^db\}
3150  if(pth(ix^d)<small_pressure) then
3151  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3152  " encountered when call mhd_get_pthermal"
3153  write(*,*) "Iteration: ", it, " Time: ", global_time
3154  write(*,*) "Location: ", x(ix^d,:)
3155  write(*,*) "Cell number: ", ix^d
3156  do iw=1,nw
3157  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3158  end do
3159  ! use erroneous arithmetic operation to crash the run
3160  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3161  write(*,*) "Saving status at the previous time step"
3162  crash=.true.
3163  end if
3164  {enddo^d&\}
3165  end if
3166 
3167  end subroutine mhd_get_pthermal_eint
3168 
3169  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
3170  subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
3173 
3174  integer, intent(in) :: ixI^L, ixO^L
3175  double precision, intent(in) :: w(ixI^S,nw)
3176  double precision, intent(in) :: x(ixI^S,1:ndim)
3177  double precision, intent(out):: pth(ixI^S)
3178  integer :: iw, ix^D
3179 
3180  pth(ixo^s)=gamma_1*(w(ixo^s,e_)&
3181  - mhd_kin_en(w,ixi^l,ixo^l)&
3182  - mhd_mag_en(w,ixi^l,ixo^l))
3183 
3184  if(has_equi_pe0) then
3185  pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe0_,b0i)
3186  end if
3187 
3188  if (fix_small_values) then
3189  {do ix^db= ixo^lim^db\}
3190  if(pth(ix^d)<small_pressure) then
3191  pth(ix^d)=small_pressure
3192  end if
3193  {enddo^d&\}
3194  else if (check_small_values) then
3195  {do ix^db= ixo^lim^db\}
3196  if(pth(ix^d)<small_pressure) then
3197  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3198  " encountered when call mhd_get_pthermal"
3199  write(*,*) "Iteration: ", it, " Time: ", global_time
3200  write(*,*) "Location: ", x(ix^d,:)
3201  write(*,*) "Cell number: ", ix^d
3202  do iw=1,nw
3203  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3204  end do
3205  ! use erroneous arithmetic operation to crash the run
3206  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3207  write(*,*) "Saving status at the previous time step"
3208  crash=.true.
3209  end if
3210  {enddo^d&\}
3211  end if
3212 
3213  end subroutine mhd_get_pthermal_origin
3214 
3215  !> Calculate thermal pressure
3216  subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
3219 
3220  integer, intent(in) :: ixI^L, ixO^L
3221  double precision, intent(in) :: w(ixI^S,nw)
3222  double precision, intent(in) :: x(ixI^S,1:ndim)
3223  double precision, intent(out):: pth(ixI^S)
3224  integer :: iw, ix^D
3225 
3226  double precision :: wprim(ixI^S,nw)
3227 
3228  wprim=w
3229  call mhd_to_primitive_semirelati(ixi^l,ixo^l,wprim,x)
3230  pth(ixo^s)=wprim(ixo^s,p_)
3231 
3232  if (.not.fix_small_values .and. check_small_values) then
3233  {do ix^db= ixo^lim^db\}
3234  if(pth(ix^d)<small_pressure) then
3235  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3236  " encountered when call mhd_get_pthermal_semirelati"
3237  write(*,*) "Iteration: ", it, " Time: ", global_time
3238  write(*,*) "Location: ", x(ix^d,:)
3239  write(*,*) "Cell number: ", ix^d
3240  do iw=1,nw
3241  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3242  end do
3243  ! use erroneous arithmetic operation to crash the run
3244  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3245  write(*,*) "Saving status at the previous time step"
3246  crash=.true.
3247  end if
3248  {enddo^d&\}
3249  end if
3250 
3251  end subroutine mhd_get_pthermal_semirelati
3252 
3253  !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
3254  subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
3257 
3258  integer, intent(in) :: ixI^L, ixO^L
3259  double precision, intent(in) :: w(ixI^S,nw)
3260  double precision, intent(in) :: x(ixI^S,1:ndim)
3261  double precision, intent(out):: pth(ixI^S)
3262  integer :: iw, ix^D
3263 
3264  pth(ixo^s)=gamma_1*(w(ixo^s,e_)-mhd_kin_en(w,ixi^l,ixo^l))
3265 
3266  if (fix_small_values) then
3267  {do ix^db= ixo^lim^db\}
3268  if(pth(ix^d)<small_pressure) then
3269  pth(ix^d)=small_pressure
3270  end if
3271  {enddo^d&\}
3272  else if (check_small_values) then
3273  {do ix^db= ixo^lim^db\}
3274  if(pth(ix^d)<small_pressure) then
3275  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3276  " encountered when call mhd_get_pthermal_hde"
3277  write(*,*) "Iteration: ", it, " Time: ", global_time
3278  write(*,*) "Location: ", x(ix^d,:)
3279  write(*,*) "Cell number: ", ix^d
3280  do iw=1,nw
3281  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3282  end do
3283  ! use erroneous arithmetic operation to crash the run
3284  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3285  write(*,*) "Saving status at the previous time step"
3286  crash=.true.
3287  end if
3288  {enddo^d&\}
3289  end if
3290 
3291  end subroutine mhd_get_pthermal_hde
3292 
3293  !> copy temperature from stored Te variable
3294  subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
3296  integer, intent(in) :: ixI^L, ixO^L
3297  double precision, intent(in) :: w(ixI^S, 1:nw)
3298  double precision, intent(in) :: x(ixI^S, 1:ndim)
3299  double precision, intent(out):: res(ixI^S)
3300  res(ixo^s) = w(ixo^s, te_)
3301  end subroutine mhd_get_temperature_from_te
3302 
3303  !> Calculate temperature=p/rho when in e_ the internal energy is stored
3304  subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
3306  integer, intent(in) :: ixI^L, ixO^L
3307  double precision, intent(in) :: w(ixI^S, 1:nw)
3308  double precision, intent(in) :: x(ixI^S, 1:ndim)
3309  double precision, intent(out):: res(ixI^S)
3310 
3311  double precision :: R(ixI^S)
3312 
3313  call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3314  res(ixo^s) = gamma_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
3315  end subroutine mhd_get_temperature_from_eint
3316 
3317  !> Calculate temperature=p/rho when in e_ the total energy is stored
3318  subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
3320  integer, intent(in) :: ixI^L, ixO^L
3321  double precision, intent(in) :: w(ixI^S, 1:nw)
3322  double precision, intent(in) :: x(ixI^S, 1:ndim)
3323  double precision, intent(out):: res(ixI^S)
3324 
3325  double precision :: R(ixI^S)
3326 
3327  call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3328  call mhd_get_pthermal(w,x,ixi^l,ixo^l,res)
3329  res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
3330 
3331  end subroutine mhd_get_temperature_from_etot
3332 
3333  subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
3335  integer, intent(in) :: ixI^L, ixO^L
3336  double precision, intent(in) :: w(ixI^S, 1:nw)
3337  double precision, intent(in) :: x(ixI^S, 1:ndim)
3338  double precision, intent(out):: res(ixI^S)
3339 
3340  double precision :: R(ixI^S)
3341 
3342  call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3343  call mhd_get_pthermal(w,x,ixi^l,ixo^l,res)
3344  res(ixo^s)=res(ixo^s)/(r(ixo^s)*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i)))
3345 
3347 
3348  subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
3350  integer, intent(in) :: ixI^L, ixO^L
3351  double precision, intent(in) :: w(ixI^S, 1:nw)
3352  double precision, intent(in) :: x(ixI^S, 1:ndim)
3353  double precision, intent(out):: res(ixI^S)
3354 
3355  double precision :: R(ixI^S)
3356 
3357  call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3358  res(ixo^s) = (gamma_1 * w(ixo^s, e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)) /&
3359  ((w(ixo^s,rho_) +block%equi_vars(ixo^s,equi_rho0_,b0i))*r(ixo^s))
3360 
3362 
3363  subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
3365  integer, intent(in) :: ixI^L, ixO^L
3366  double precision, intent(in) :: w(ixI^S, 1:nw)
3367  double precision, intent(in) :: x(ixI^S, 1:ndim)
3368  double precision, intent(out):: res(ixI^S)
3369 
3370  double precision :: R(ixI^S)
3371 
3372  call mhd_get_rfactor(w,x,ixi^l,ixo^l,r)
3373  res(ixo^s)= block%equi_vars(ixo^s,equi_pe0_,b0i)/(block%equi_vars(ixo^s,equi_rho0_,b0i)*r(ixo^s))
3374 
3375  end subroutine mhd_get_temperature_equi
3376 
3377  subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
3379  integer, intent(in) :: ixI^L, ixO^L
3380  double precision, intent(in) :: w(ixI^S, 1:nw)
3381  double precision, intent(in) :: x(ixI^S, 1:ndim)
3382  double precision, intent(out):: res(ixI^S)
3383  res(ixo^s) = block%equi_vars(ixo^s,equi_rho0_,b0i)
3384  end subroutine mhd_get_rho_equi
3385 
3386  subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
3388  integer, intent(in) :: ixI^L, ixO^L
3389  double precision, intent(in) :: w(ixI^S, 1:nw)
3390  double precision, intent(in) :: x(ixI^S, 1:ndim)
3391  double precision, intent(out):: res(ixI^S)
3392  res(ixo^s) = block%equi_vars(ixo^s,equi_pe0_,b0i)
3393  end subroutine mhd_get_pe_equi
3394 
3395  !> Calculate the square of the thermal sound speed csound2 within ixO^L.
3396  !> csound2=gamma*p/rho
3397  subroutine mhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
3399  integer, intent(in) :: ixi^l, ixo^l
3400  double precision, intent(in) :: w(ixi^s,nw)
3401  double precision, intent(in) :: x(ixi^s,1:ndim)
3402  double precision, intent(out) :: csound2(ixi^s)
3403  double precision :: rho(ixi^s)
3404 
3405  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
3406  if(mhd_energy) then
3407  call mhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
3408  csound2(ixo^s)=mhd_gamma*csound2(ixo^s)/rho(ixo^s)
3409  else
3410  csound2(ixo^s)=mhd_gamma*mhd_adiab*rho(ixo^s)**gamma_1
3411  end if
3412  end subroutine mhd_get_csound2
3413 
3414  !> Calculate total pressure within ixO^L including magnetic pressure
3415  subroutine mhd_get_p_total(w,x,ixI^L,ixO^L,p)
3417 
3418  integer, intent(in) :: ixI^L, ixO^L
3419  double precision, intent(in) :: w(ixI^S,nw)
3420  double precision, intent(in) :: x(ixI^S,1:ndim)
3421  double precision, intent(out) :: p(ixI^S)
3422 
3423  call mhd_get_pthermal(w,x,ixi^l,ixo^l,p)
3424 
3425  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3426 
3427  end subroutine mhd_get_p_total
3428 
3429  !> Calculate fluxes within ixO^L without any splitting
3430  subroutine mhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3432  use mod_geometry
3433 
3434  integer, intent(in) :: ixI^L, ixO^L, idim
3435  ! conservative w
3436  double precision, intent(in) :: wC(ixI^S,nw)
3437  ! primitive w
3438  double precision, intent(in) :: w(ixI^S,nw)
3439  double precision, intent(in) :: x(ixI^S,1:ndim)
3440  double precision,intent(out) :: f(ixI^S,nwflux)
3441 
3442  double precision :: ptotal(ixO^S)
3443  double precision :: tmp(ixI^S)
3444  double precision :: vHall(ixI^S,1:ndir)
3445  integer :: idirmin, iw, idir, jdir, kdir
3446  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3447  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3448 
3449  ! Get flux of density
3450  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3451 
3452  if(mhd_energy) then
3453  ptotal(ixo^s)=w(ixo^s,p_)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3454  else
3455  ptotal(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3456  end if
3457 
3458  if (mhd_hall) then
3459  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3460  end if
3461 
3462  ! Get flux of tracer
3463  do iw=1,mhd_n_tracer
3464  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3465  end do
3466 
3467  ! Get flux of momentum
3468  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3469  do idir=1,ndir
3470  if(idim==idir) then
3471  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3472  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3473  else
3474  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3475  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3476  end if
3477  end do
3478 
3479  ! Get flux of energy
3480  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
3481  if(mhd_energy) then
3482  if (mhd_internal_e) then
3483  f(ixo^s,e_)=w(ixo^s,mom(idim))*wc(ixo^s,e_)
3484  else
3485  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
3486  -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
3487  if(mhd_hall) then
3488  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
3489  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
3490  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3491  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3492  end if
3493  end if
3494  end if
3495 
3496  ! compute flux of magnetic field
3497  ! f_i[b_k]=v_i*b_k-v_k*b_i
3498  do idir=1,ndir
3499  if (idim==idir) then
3500  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3501  if (mhd_glm) then
3502  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3503  else
3504  f(ixo^s,mag(idir))=zero
3505  end if
3506  else
3507  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
3508  if (mhd_hall) then
3509  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
3510  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3511  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3512  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3513  end if
3514  end if
3515  end do
3516 
3517  if (mhd_glm) then
3518  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3519  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3520  end if
3521 
3522  ! Contributions of ambipolar term in explicit scheme
3523  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3524  ! ambipolar electric field
3525  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3526  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3527  allocate(jambi(ixi^s,1:3))
3528  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3529  allocate(btot(ixo^s,1:3))
3530  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3531  allocate(tmp2(ixo^s),tmp3(ixo^s))
3532  !tmp2 = Btot^2
3533  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3534  !tmp3 = J_ambi dot Btot
3535  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3536 
3537  select case(idim)
3538  case(1)
3539  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3540  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3541  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3542  case(2)
3543  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3544  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3545  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3546  case(3)
3547  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3548  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3549  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3550  endselect
3551 
3552  if(mhd_energy .and. .not. mhd_internal_e) then
3553  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3554  endif
3555 
3556  deallocate(jambi,btot,tmp2,tmp3)
3557  endif
3558 
3559  end subroutine mhd_get_flux
3560 
3561  !> Calculate fluxes within ixO^L without any splitting
3562  subroutine mhd_get_flux_hde(wC,w,x,ixI^L,ixO^L,idim,f)
3564  use mod_geometry
3565 
3566  integer, intent(in) :: ixI^L, ixO^L, idim
3567  ! conservative w
3568  double precision, intent(in) :: wC(ixI^S,nw)
3569  ! primitive w
3570  double precision, intent(in) :: w(ixI^S,nw)
3571  double precision, intent(in) :: x(ixI^S,1:ndim)
3572  double precision,intent(out) :: f(ixI^S,nwflux)
3573 
3574  double precision :: pgas(ixO^S), ptotal(ixO^S)
3575  double precision :: tmp(ixI^S)
3576  integer :: idirmin, iw, idir, jdir, kdir
3577  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3578  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3579 
3580  ! Get flux of density
3581  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3582  ! pgas is time dependent only
3583  if(mhd_energy) then
3584  pgas(ixo^s)=w(ixo^s,p_)
3585  else
3586  pgas(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3587  end if
3588 
3589  ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3590 
3591  ! Get flux of tracer
3592  do iw=1,mhd_n_tracer
3593  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3594  end do
3595 
3596  ! Get flux of momentum
3597  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3598  do idir=1,ndir
3599  if(idim==idir) then
3600  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3601  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3602  else
3603  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3604  w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3605  end if
3606  end do
3607 
3608  ! Get flux of energy
3609  if(mhd_energy) then
3610  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+pgas(ixo^s))
3611  end if
3612 
3613  ! compute flux of magnetic field
3614  ! f_i[b_k]=v_i*b_k-v_k*b_i
3615  do idir=1,ndir
3616  if (idim==idir) then
3617  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3618  if (mhd_glm) then
3619  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3620  else
3621  f(ixo^s,mag(idir))=zero
3622  end if
3623  else
3624  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
3625  end if
3626  end do
3627 
3628  if (mhd_glm) then
3629  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3630  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3631  end if
3632 
3633  ! Contributions of ambipolar term in explicit scheme
3634  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3635  ! ambipolar electric field
3636  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3637  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3638  allocate(jambi(ixi^s,1:3))
3639  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3640  allocate(btot(ixo^s,1:3))
3641  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3642  allocate(tmp2(ixo^s),tmp3(ixo^s))
3643  !tmp2 = Btot^2
3644  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3645  !tmp3 = J_ambi dot Btot
3646  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3647 
3648  select case(idim)
3649  case(1)
3650  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3651  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3652  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3653  case(2)
3654  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3655  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3656  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3657  case(3)
3658  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3659  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3660  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3661  endselect
3662 
3663  if(mhd_energy) then
3664  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3665  endif
3666 
3667  deallocate(jambi,btot,tmp2,tmp3)
3668  endif
3669 
3670  end subroutine mhd_get_flux_hde
3671 
3672  !> Calculate fluxes within ixO^L with possible splitting
3673  subroutine mhd_get_flux_split(wC,w,x,ixI^L,ixO^L,idim,f)
3675  use mod_geometry
3676 
3677  integer, intent(in) :: ixI^L, ixO^L, idim
3678  ! conservative w
3679  double precision, intent(in) :: wC(ixI^S,nw)
3680  ! primitive w
3681  double precision, intent(in) :: w(ixI^S,nw)
3682  double precision, intent(in) :: x(ixI^S,1:ndim)
3683  double precision,intent(out) :: f(ixI^S,nwflux)
3684 
3685  double precision :: pgas(ixO^S), ptotal(ixO^S), B(ixO^S,1:ndir)
3686  double precision :: tmp(ixI^S)
3687  double precision :: vHall(ixI^S,1:ndir)
3688  integer :: idirmin, iw, idir, jdir, kdir
3689  double precision, allocatable, dimension(:^D&,:) :: Jambi, btot
3690  double precision, allocatable, dimension(:^D&) :: tmp2, tmp3
3691  double precision :: tmp4(ixO^S)
3692 
3693 
3694  call mhd_get_rho(w,x,ixi^l,ixo^l,tmp)
3695  ! Get flux of density
3696  f(ixo^s,rho_)=w(ixo^s,mom(idim))*tmp(ixo^s)
3697  ! pgas is time dependent only
3698  if(mhd_energy) then
3699  pgas(ixo^s)=w(ixo^s,p_)
3700  else
3701  pgas(ixo^s)=mhd_adiab*tmp(ixo^s)**mhd_gamma
3702  if(has_equi_pe0) then
3703  pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe0_,b0i)
3704  endif
3705  end if
3706 
3707  if (mhd_hall) then
3708  call mhd_getv_hall(w,x,ixi^l,ixo^l,vhall)
3709  end if
3710 
3711  if(b0field) then
3712  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,idim)
3713  pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,idim),dim=ndim+1)
3714  else
3715  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3716  end if
3717 
3718  ptotal(ixo^s)=pgas(ixo^s)+0.5d0*sum(w(ixo^s,mag(:))**2,dim=ndim+1)
3719 
3720  ! Get flux of tracer
3721  do iw=1,mhd_n_tracer
3722  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3723  end do
3724 
3725  ! Get flux of momentum
3726  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3727  if(b0field) then
3728  do idir=1,ndir
3729  if(idim==idir) then
3730  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3731  w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3732  block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3733  else
3734  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3735  w(ixo^s,mag(idir))*b(ixo^s,idim)-&
3736  block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3737  end if
3738  end do
3739  else
3740  do idir=1,ndir
3741  if(idim==idir) then
3742  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))+ptotal(ixo^s)-&
3743  w(ixo^s,mag(idir))*b(ixo^s,idim)
3744  else
3745  f(ixo^s,mom(idir))=wc(ixo^s,mom(idir))*w(ixo^s,mom(idim))-&
3746  w(ixo^s,mag(idir))*b(ixo^s,idim)
3747  end if
3748  end do
3749  end if
3750 
3751  ! Get flux of energy
3752  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
3753  if(mhd_energy) then
3754  if (mhd_internal_e) then
3755  f(ixo^s,e_)=w(ixo^s,mom(idim))*wc(ixo^s,e_)
3756  else
3757  f(ixo^s,e_)=w(ixo^s,mom(idim))*(wc(ixo^s,e_)+ptotal(ixo^s))&
3758  -b(ixo^s,idim)*sum(w(ixo^s,mag(:))*w(ixo^s,mom(:)),dim=ndim+1)
3759 
3760  if (mhd_hall) then
3761  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
3762  if (mhd_etah>zero) then
3763  f(ixo^s,e_) = f(ixo^s,e_) + vhall(ixo^s,idim) * &
3764  sum(w(ixo^s, mag(:))*b(ixo^s,:),dim=ndim+1) &
3765  - b(ixo^s,idim) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3766  end if
3767  end if
3768  end if
3769  if(has_equi_pe0) then
3770  f(ixo^s,e_)= f(ixo^s,e_) &
3771  + w(ixo^s,mom(idim)) * block%equi_vars(ixo^s,equi_pe0_,idim) * inv_gamma_1
3772  end if
3773  end if
3774 
3775  ! compute flux of magnetic field
3776  ! f_i[b_k]=v_i*b_k-v_k*b_i
3777  do idir=1,ndir
3778  if (idim==idir) then
3779  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3780  if (mhd_glm) then
3781  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3782  else
3783  f(ixo^s,mag(idir))=zero
3784  end if
3785  else
3786  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,mom(idir))
3787 
3788  if (mhd_hall) then
3789  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
3790  if (mhd_etah>zero) then
3791  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3792  - vhall(ixo^s,idir)*b(ixo^s,idim) &
3793  + vhall(ixo^s,idim)*b(ixo^s,idir)
3794  end if
3795  end if
3796 
3797  end if
3798  end do
3799 
3800  if (mhd_glm) then
3801  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3802  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3803  end if
3804 
3805  ! Contributions of ambipolar term in explicit scheme
3806  if(mhd_ambipolar_exp.and. .not.stagger_grid) then
3807  ! ambipolar electric field
3808  ! E_ambi=-eta_ambi*JxBxB=-JaxBxB=B^2*Ja-(Ja dot B)*B
3809  !Ja=eta_ambi*J=J * mhd_eta_ambi/rho**2
3810  allocate(jambi(ixi^s,1:3))
3811  call mhd_get_jambi(w,x,ixi^l,ixo^l,jambi)
3812  allocate(btot(ixo^s,1:3))
3813  if(b0field) then
3814  do idir=1,3
3815  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,idim)
3816  enddo
3817  else
3818  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
3819  endif
3820  allocate(tmp2(ixo^s),tmp3(ixo^s))
3821  !tmp2 = Btot^2
3822  tmp2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1)
3823  !tmp3 = J_ambi dot Btot
3824  tmp3(ixo^s) = sum(jambi(ixo^s,:)*btot(ixo^s,:),dim=ndim+1)
3825 
3826  select case(idim)
3827  case(1)
3828  tmp(ixo^s)=w(ixo^s,mag(3)) *jambi(ixo^s,2) - w(ixo^s,mag(2)) * jambi(ixo^s,3)
3829  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(2)) * btot(ixo^s,3) - w(ixo^s,mag(3)) * btot(ixo^s,2)
3830  f(ixo^s,mag(2))= f(ixo^s,mag(2)) - tmp2(ixo^s) * jambi(ixo^s,3) + tmp3(ixo^s) * btot(ixo^s,3)
3831  f(ixo^s,mag(3))= f(ixo^s,mag(3)) + tmp2(ixo^s) * jambi(ixo^s,2) - tmp3(ixo^s) * btot(ixo^s,2)
3832  case(2)
3833  tmp(ixo^s)=w(ixo^s,mag(1)) *jambi(ixo^s,3) - w(ixo^s,mag(3)) * jambi(ixo^s,1)
3834  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(3)) * btot(ixo^s,1) - w(ixo^s,mag(1)) * btot(ixo^s,3)
3835  f(ixo^s,mag(1))= f(ixo^s,mag(1)) + tmp2(ixo^s) * jambi(ixo^s,3) - tmp3(ixo^s) * btot(ixo^s,3)
3836  f(ixo^s,mag(3))= f(ixo^s,mag(3)) - tmp2(ixo^s) * jambi(ixo^s,1) + tmp3(ixo^s) * btot(ixo^s,1)
3837  case(3)
3838  tmp(ixo^s)=w(ixo^s,mag(2)) *jambi(ixo^s,1) - w(ixo^s,mag(1)) * jambi(ixo^s,2)
3839  if(b0field) tmp4(ixo^s) = w(ixo^s,mag(1)) * btot(ixo^s,2) - w(ixo^s,mag(2)) * btot(ixo^s,1)
3840  f(ixo^s,mag(1))= f(ixo^s,mag(1)) - tmp2(ixo^s) * jambi(ixo^s,2) + tmp3(ixo^s) * btot(ixo^s,2)
3841  f(ixo^s,mag(2))= f(ixo^s,mag(2)) + tmp2(ixo^s) * jambi(ixo^s,1) - tmp3(ixo^s) * btot(ixo^s,1)
3842  endselect
3843 
3844  if(mhd_energy .and. .not. mhd_internal_e) then
3845  f(ixo^s,e_) = f(ixo^s,e_) + tmp2(ixo^s) * tmp(ixo^s)
3846  if(b0field) f(ixo^s,e_) = f(ixo^s,e_) + tmp3(ixo^s) * tmp4(ixo^s)
3847  endif
3848 
3849  deallocate(jambi,btot,tmp2,tmp3)
3850  endif
3851 
3852  end subroutine mhd_get_flux_split
3853 
3854  !> Calculate semirelativistic fluxes within ixO^L without any splitting
3855  subroutine mhd_get_flux_semirelati(wC,w,x,ixI^L,ixO^L,idim,f)
3857  use mod_geometry
3858 
3859  integer, intent(in) :: ixI^L, ixO^L, idim
3860  ! conservative w
3861  double precision, intent(in) :: wC(ixI^S,nw)
3862  ! primitive w
3863  double precision, intent(in) :: w(ixI^S,nw)
3864  double precision, intent(in) :: x(ixI^S,1:ndim)
3865  double precision,intent(out) :: f(ixI^S,nwflux)
3866 
3867  double precision :: pgas(ixO^S)
3868  double precision :: SA(ixO^S), E(ixO^S,1:ndir), B(ixO^S,1:ndir)
3869  integer :: idirmin, iw, idir, jdir, kdir
3870 
3871  ! gas thermal pressure
3872  if(mhd_energy) then
3873  pgas(ixo^s)=w(ixo^s,p_)
3874  else
3875  pgas(ixo^s)=mhd_adiab*w(ixo^s,rho_)**mhd_gamma
3876  end if
3877 
3878  ! Get flux of density
3879  f(ixo^s,rho_)=w(ixo^s,mom(idim))*w(ixo^s,rho_)
3880 
3881  ! Get flux of tracer
3882  do iw=1,mhd_n_tracer
3883  f(ixo^s,tracer(iw))=w(ixo^s,mom(idim))*w(ixo^s,tracer(iw))
3884  end do
3885  ! E=-uxB=Bxu
3886  if(b0field) then
3887  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))+block%B0(ixo^s,1:ndir,idim)
3888  pgas(ixo^s)=pgas(ixo^s)+sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,idim),dim=ndim+1)
3889  else
3890  b(ixo^s,1:ndir)=w(ixo^s,mag(1:ndir))
3891  end if
3892  e=0.d0
3893  do idir=1,ndir; do jdir=1,ndir; do kdir=1,ndir
3894  if(lvc(idir,jdir,kdir)==1)then
3895  e(ixo^s,idir)=e(ixo^s,idir)+b(ixo^s,jdir)*w(ixo^s,mom(kdir))
3896  else if(lvc(idir,jdir,kdir)==-1)then
3897  e(ixo^s,idir)=e(ixo^s,idir)-b(ixo^s,jdir)*w(ixo^s,mom(kdir))
3898  end if
3899  end do; end do; end do
3900 
3901  pgas(ixo^s)=pgas(ixo^s)+half*(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+&
3902  sum(e(ixo^s,:)**2,dim=ndim+1)*inv_squared_c)
3903 
3904  ! Get flux of momentum
3905  if(b0field) then
3906  do idir=1,ndir
3907  if(idim==idir) then
3908  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))+pgas&
3909  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3910  -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3911  else
3912  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))&
3913  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c&
3914  -block%B0(ixo^s,idir,idim)*w(ixo^s,mag(idim))
3915  end if
3916  end do
3917  else
3918  do idir=1,ndir
3919  if(idim==idir) then
3920  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))+pgas&
3921  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3922  else
3923  f(ixo^s,mom(idir))=w(ixo^s,rho_)*w(ixo^s,mom(idir))*w(ixo^s,mom(idim))&
3924  -w(ixo^s,mag(idir))*b(ixo^s,idim)-e(ixo^s,idir)*e(ixo^s,idim)*inv_squared_c
3925  end if
3926  end do
3927  end if
3928 
3929  ! Get flux of total energy
3930  if(mhd_internal_e) then
3931  ! Get flux of internal energy
3932  f(ixo^s,e_)=w(ixo^s,mom(idim))*wc(ixo^s,e_)
3933  else if(mhd_energy) then
3934  sa=0.d0
3935  do jdir=1,ndir; do kdir=1,ndir
3936  if(lvc(idim,jdir,kdir)==1)then
3937  sa(ixo^s)=sa(ixo^s)+e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3938  else if(lvc(idim,jdir,kdir)==-1) then
3939  sa(ixo^s)=sa(ixo^s)-e(ixo^s,jdir)*w(ixo^s,mag(kdir))
3940  end if
3941  end do; end do
3942  f(ixo^s,e_)=w(ixo^s,mom(idim))*(half*w(ixo^s,rho_)*sum(w(ixo^s,mom(:))**2,dim=ndim+1)+&
3943  mhd_gamma*pgas*inv_gamma_1)+sa(ixo^s)
3944  end if
3945 
3946  ! compute flux of magnetic field
3947  ! f_i[b_k]=v_i*b_k-v_k*b_i
3948  do idir=1,ndir
3949  if (idim==idir) then
3950  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3951  if (mhd_glm) then
3952  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3953  else
3954  f(ixo^s,mag(idir))=zero
3955  end if
3956  else
3957  f(ixo^s,mag(idir))=w(ixo^s,mom(idim))*b(ixo^s,idir)-b(ixo^s,idim)*w(ixo^s,mom(idir))
3958  end if
3959  end do
3960 
3961  if (mhd_glm) then
3962  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3963  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3964  end if
3965 
3966  end subroutine mhd_get_flux_semirelati
3967 
3968  !> Source terms J.E in internal energy.
3969  !> For the ambipolar term E = ambiCoef * JxBxB=ambiCoef * B^2(-J_perpB)
3970  !=> the source term J.E = ambiCoef * B^2 * J_perpB^2 = ambiCoef * JxBxB^2/B^2
3971  !> ambiCoef is calculated as mhd_ambi_coef/rho^2, see also the subroutine mhd_get_Jambi
3972  subroutine add_source_ambipolar_internal_energy(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3974  integer, intent(in) :: ixI^L, ixO^L,ie
3975  double precision, intent(in) :: qdt
3976  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3977  double precision, intent(inout) :: w(ixI^S,1:nw)
3978  double precision :: tmp(ixI^S)
3979  double precision :: jxbxb(ixI^S,1:3)
3980 
3981  call mhd_get_jxbxb(wct,x,ixi^l,ixo^l,jxbxb)
3982  tmp(ixo^s) = sum(jxbxb(ixo^s,1:3)**2,dim=ndim+1) / mhd_mag_en_all(wct, ixi^l, ixo^l)
3983  call multiplyambicoef(ixi^l,ixo^l,tmp,wct,x)
3984  w(ixo^s,ie)=w(ixo^s,ie)+qdt * tmp
3985 
3987 
3988  subroutine mhd_get_jxbxb(w,x,ixI^L,ixO^L,res)
3990 
3991  integer, intent(in) :: ixI^L, ixO^L
3992  double precision, intent(in) :: w(ixI^S,nw)
3993  double precision, intent(in) :: x(ixI^S,1:ndim)
3994  double precision, intent(out) :: res(:^D&,:)
3995 
3996  double precision :: btot(ixI^S,1:3)
3997  integer :: idir, idirmin
3998  double precision :: current(ixI^S,7-2*ndir:3)
3999  double precision :: tmp(ixI^S),b2(ixI^S)
4000 
4001  res=0.d0
4002  ! Calculate current density and idirmin
4003  call get_current(w,ixi^l,ixo^l,idirmin,current)
4004  !!!here we know that current has nonzero values only for components in the range idirmin, 3
4005 
4006  if(b0field) then
4007  do idir=1,3
4008  btot(ixo^s, idir) = w(ixo^s,mag(idir)) + block%B0(ixo^s,idir,b0i)
4009  enddo
4010  else
4011  btot(ixo^s,1:3) = w(ixo^s,mag(1:3))
4012  endif
4013 
4014  tmp(ixo^s) = sum(current(ixo^s,idirmin:3)*btot(ixo^s,idirmin:3),dim=ndim+1) !J.B
4015  b2(ixo^s) = sum(btot(ixo^s,1:3)**2,dim=ndim+1) !B^2
4016  do idir=1,idirmin-1
4017  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s)
4018  enddo
4019  do idir=idirmin,3
4020  res(ixo^s,idir) = btot(ixo^s,idir) * tmp(ixo^s) - current(ixo^s,idir) * b2(ixo^s)
4021  enddo
4022  end subroutine mhd_get_jxbxb
4023 
4024  !> Sets the sources for the ambipolar
4025  !> this is used for the STS method
4026  ! The sources are added directly (instead of fluxes as in the explicit)
4027  !> at the corresponding indices
4028  !> store_flux_var is explicitly called for each of the fluxes one by one
4029  subroutine sts_set_source_ambipolar(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
4031  use mod_fix_conserve
4032 
4033  integer, intent(in) :: ixI^L, ixO^L,igrid,nflux
4034  double precision, intent(in) :: x(ixI^S,1:ndim)
4035  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
4036  double precision, intent(in) :: my_dt
4037  logical, intent(in) :: fix_conserve_at_step
4038 
4039  double precision, dimension(ixI^S,1:3) :: tmp,ff
4040  double precision :: fluxall(ixI^S,1:nflux,1:ndim)
4041  double precision :: fE(ixI^S,sdim:3)
4042  double precision :: btot(ixI^S,1:3),tmp2(ixI^S)
4043  integer :: i, ixA^L, ie_
4044 
4045  ixa^l=ixo^l^ladd1;
4046 
4047  fluxall=zero
4048 
4049  call mhd_get_jxbxb(w,x,ixi^l,ixa^l,tmp)
4050 
4051  !set electric field in tmp: E=nuA * jxbxb, where nuA=-etaA/rho^2
4052  do i=1,3
4053  !tmp(ixA^S,i) = -(mhd_eta_ambi/w(ixA^S, rho_)**2) * tmp(ixA^S,i)
4054  call multiplyambicoef(ixi^l,ixa^l,tmp(ixi^s,i),w,x)
4055  enddo
4056 
4057  if(mhd_energy .and. .not.mhd_internal_e) then
4058  !btot should be only mag. pert.
4059  btot(ixa^s,1:3)=0.d0
4060  !if(B0field) then
4061  ! do i=1,ndir
4062  ! btot(ixA^S, i) = w(ixA^S,mag(i)) + block%B0(ixA^S,i,0)
4063  ! enddo
4064  !else
4065  btot(ixa^s,1:ndir) = w(ixa^s,mag(1:ndir))
4066  !endif
4067  call cross_product(ixi^l,ixa^l,tmp,btot,ff)
4068  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
4069  if(fix_conserve_at_step) fluxall(ixi^s,1,1:ndim)=ff(ixi^s,1:ndim)
4070  !- sign comes from the fact that the flux divergence is a source now
4071  wres(ixo^s,e_)=-tmp2(ixo^s)
4072  endif
4073 
4074  if(stagger_grid) then
4075  if(ndir>ndim) then
4076  !!!Bz
4077  ff(ixa^s,1) = tmp(ixa^s,2)
4078  ff(ixa^s,2) = -tmp(ixa^s,1)
4079  ff(ixa^s,3) = 0.d0
4080  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
4081  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
4082  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
4083  end if
4084  fe=0.d0
4085  call update_faces_ambipolar(ixi^l,ixo^l,w,x,tmp,fe,btot)
4086  ixamax^d=ixomax^d;
4087  ixamin^d=ixomin^d-1;
4088  wres(ixa^s,mag(1:ndim))=-btot(ixa^s,1:ndim)
4089  else
4090  !write curl(ele) as the divergence
4091  !m1={0,ele[[3]],-ele[[2]]}
4092  !m2={-ele[[3]],0,ele[[1]]}
4093  !m3={ele[[2]],-ele[[1]],0}
4094 
4095  !!!Bx
4096  ff(ixa^s,1) = 0.d0
4097  ff(ixa^s,2) = tmp(ixa^s,3)
4098  ff(ixa^s,3) = -tmp(ixa^s,2)
4099  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
4100  if(fix_conserve_at_step) fluxall(ixi^s,2,1:ndim)=ff(ixi^s,1:ndim)
4101  !flux divergence is a source now
4102  wres(ixo^s,mag(1))=-tmp2(ixo^s)
4103  !!!By
4104  ff(ixa^s,1) = -tmp(ixa^s,3)
4105  ff(ixa^s,2) = 0.d0
4106  ff(ixa^s,3) = tmp(ixa^s,1)
4107  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
4108  if(fix_conserve_at_step) fluxall(ixi^s,3,1:ndim)=ff(ixi^s,1:ndim)
4109  wres(ixo^s,mag(2))=-tmp2(ixo^s)
4110 
4111  if(ndir==3) then
4112  !!!Bz
4113  ff(ixa^s,1) = tmp(ixa^s,2)
4114  ff(ixa^s,2) = -tmp(ixa^s,1)
4115  ff(ixa^s,3) = 0.d0
4116  call get_flux_on_cell_face(ixi^l,ixo^l,ff,tmp2)
4117  if(fix_conserve_at_step) fluxall(ixi^s,1+ndir,1:ndim)=ff(ixi^s,1:ndim)
4118  wres(ixo^s,mag(ndir))=-tmp2(ixo^s)
4119  end if
4120 
4121  end if
4122 
4123  if(fix_conserve_at_step) then
4124  fluxall=my_dt*fluxall
4125  call store_flux(igrid,fluxall,1,ndim,nflux)
4126  if(stagger_grid) then
4127  call store_edge(igrid,ixi^l,my_dt*fe,1,ndim)
4128  end if
4129  end if
4130 
4131  end subroutine sts_set_source_ambipolar
4132 
4133  !> get ambipolar electric field and the integrals around cell faces
4134  subroutine update_faces_ambipolar(ixI^L,ixO^L,w,x,ECC,fE,circ)
4136 
4137  integer, intent(in) :: ixI^L, ixO^L
4138  double precision, intent(in) :: w(ixI^S,1:nw)
4139  double precision, intent(in) :: x(ixI^S,1:ndim)
4140  ! amibipolar electric field at cell centers
4141  double precision, intent(in) :: ECC(ixI^S,1:3)
4142  double precision, intent(out) :: fE(ixI^S,sdim:3)
4143  double precision, intent(out) :: circ(ixI^S,1:ndim)
4144 
4145  integer :: hxC^L,ixC^L,ixA^L
4146  integer :: idim1,idim2,idir,ix^D
4147 
4148  fe=zero
4149  ! calcuate ambipolar electric field on cell edges from cell centers
4150  do idir=sdim,3
4151  ixcmax^d=ixomax^d;
4152  ixcmin^d=ixomin^d+kr(idir,^d)-1;
4153  {do ix^db=0,1\}
4154  if({ ix^d==1 .and. ^d==idir | .or.}) cycle
4155  ixamin^d=ixcmin^d+ix^d;
4156  ixamax^d=ixcmax^d+ix^d;
4157  fe(ixc^s,idir)=fe(ixc^s,idir)+ecc(ixa^s,idir)
4158  {end do\}
4159  fe(ixc^s,idir)=fe(ixc^s,idir)*0.25d0*block%dsC(ixc^s,idir)
4160  end do
4161 
4162  ! Calculate circulation on each face to get value of line integral of
4163  ! electric field in the positive idir direction.
4164  ixcmax^d=ixomax^d;
4165  ixcmin^d=ixomin^d-1;
4166 
4167  circ=zero
4168 
4169  do idim1=1,ndim ! Coordinate perpendicular to face
4170  do idim2=1,ndim
4171  do idir=sdim,3 ! Direction of line integral
4172  ! Assemble indices
4173  hxc^l=ixc^l-kr(idim2,^d);
4174  ! Add line integrals in direction idir
4175  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
4176  +lvc(idim1,idim2,idir)&
4177  *(fe(ixc^s,idir)&
4178  -fe(hxc^s,idir))
4179  end do
4180  end do
4181  circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
4182  end do
4183 
4184  end subroutine update_faces_ambipolar
4185 
4186  !> use cell-center flux to get cell-face flux
4187  !> and get the source term as the divergence of the flux
4188  subroutine get_flux_on_cell_face(ixI^L,ixO^L,ff,src)
4190 
4191  integer, intent(in) :: ixI^L, ixO^L
4192  double precision, dimension(:^D&,:), intent(inout) :: ff
4193  double precision, intent(out) :: src(ixI^S)
4194 
4195  double precision :: ffc(ixI^S,1:ndim)
4196  double precision :: dxinv(ndim)
4197  integer :: idims, ix^D, ixA^L, ixB^L, ixC^L
4198 
4199  ixa^l=ixo^l^ladd1;
4200  dxinv=1.d0/dxlevel
4201  ! cell corner flux in ffc
4202  ffc=0.d0
4203  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
4204  {do ix^db=0,1\}
4205  ixbmin^d=ixcmin^d+ix^d;
4206  ixbmax^d=ixcmax^d+ix^d;
4207  ffc(ixc^s,1:ndim)=ffc(ixc^s,1:ndim)+ff(ixb^s,1:ndim)
4208  {end do\}
4209  ffc(ixc^s,1:ndim)=0.5d0**ndim*ffc(ixc^s,1:ndim)
4210  ! flux at cell face
4211  ff(ixi^s,1:ndim)=0.d0
4212  do idims=1,ndim
4213  ixb^l=ixo^l-kr(idims,^d);
4214  ixcmax^d=ixomax^d; ixcmin^d=ixbmin^d;
4215  {do ix^db=0,1 \}
4216  if({ ix^d==0 .and. ^d==idims | .or.}) then
4217  ixbmin^d=ixcmin^d-ix^d;
4218  ixbmax^d=ixcmax^d-ix^d;
4219  ff(ixc^s,idims)=ff(ixc^s,idims)+ffc(ixb^s,idims)
4220  end if
4221  {end do\}
4222  ff(ixc^s,idims)=ff(ixc^s,idims)*0.5d0**(ndim-1)
4223  end do
4224  src=0.d0
4225  if(slab_uniform) then
4226  do idims=1,ndim
4227  ff(ixa^s,idims)=dxinv(idims)*ff(ixa^s,idims)
4228  ixb^l=ixo^l-kr(idims,^d);
4229  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4230  end do
4231  else
4232  do idims=1,ndim
4233  ff(ixa^s,idims)=ff(ixa^s,idims)*block%surfaceC(ixa^s,idims)
4234  ixb^l=ixo^l-kr(idims,^d);
4235  src(ixo^s)=src(ixo^s)+ff(ixo^s,idims)-ff(ixb^s,idims)
4236  end do
4237  src(ixo^s)=src(ixo^s)/block%dvolume(ixo^s)
4238  end if
4239  end subroutine get_flux_on_cell_face
4240 
4241  !> Calculates the explicit dt for the ambipokar term
4242  !> This function is used by both explicit scheme and STS method
4243  function get_ambipolar_dt(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
4245 
4246  integer, intent(in) :: ixi^l, ixo^l
4247  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
4248  double precision, intent(in) :: w(ixi^s,1:nw)
4249  double precision :: dtnew
4250 
4251  double precision :: coef
4252  double precision :: dxarr(ndim)
4253  double precision :: tmp(ixi^s)
4254 
4255  ^d&dxarr(^d)=dx^d;
4256  tmp(ixo^s) = mhd_mag_en_all(w, ixi^l, ixo^l)
4257  call multiplyambicoef(ixi^l,ixo^l,tmp,w,x)
4258  coef = maxval(abs(tmp(ixo^s)))
4259  if(coef/=0.d0) then
4260  coef=1.d0/coef
4261  else
4262  coef=bigdouble
4263  end if
4264  if(slab_uniform) then
4265  dtnew=minval(dxarr(1:ndim))**2.0d0*coef
4266  else
4267  dtnew=minval(block%ds(ixo^s,1:ndim))**2.0d0*coef
4268  end if
4269 
4270  end function get_ambipolar_dt
4271 
4272  !> multiply res by the ambipolar coefficient
4273  !> The ambipolar coefficient is calculated as -mhd_eta_ambi/rho^2
4274  !> The user may mask its value in the user file
4275  !> by implemneting usr_mask_ambipolar subroutine
4276  subroutine multiplyambicoef(ixI^L,ixO^L,res,w,x)
4278  integer, intent(in) :: ixi^l, ixo^l
4279  double precision, intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
4280  double precision, intent(inout) :: res(ixi^s)
4281  double precision :: tmp(ixi^s)
4282  double precision :: rho(ixi^s)
4283 
4284  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
4285  tmp=0.d0
4286  tmp(ixo^s)=-mhd_eta_ambi/rho(ixo^s)**2
4287  if (associated(usr_mask_ambipolar)) then
4288  call usr_mask_ambipolar(ixi^l,ixo^l,w,x,tmp)
4289  end if
4290 
4291  res(ixo^s) = tmp(ixo^s) * res(ixo^s)
4292  end subroutine multiplyambicoef
4293 
4294  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
4295  subroutine mhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
4299  use mod_gravity, only: gravity_add_source
4300  use mod_cak_force, only: cak_add_source
4301 
4302  integer, intent(in) :: ixI^L, ixO^L
4303  double precision, intent(in) :: qdt,dtfactor
4304  double precision, intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
4305  double precision, intent(inout) :: w(ixI^S,1:nw)
4306  logical, intent(in) :: qsourcesplit
4307  logical, intent(inout) :: active
4308 
4309  !TODO local_timestep support is only added for splitting
4310  ! but not for other nonideal terms such gravity, RC, viscosity,..
4311  ! it will also only work for divbfix 'linde', which does not require
4312  ! modification as it does not use dt in the update
4313 
4314  if (.not. qsourcesplit) then
4315  if(mhd_internal_e) then
4316  ! Source for solving internal energy
4317  active = .true.
4318  call add_source_internal_e(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
4319  else
4320  if(has_equi_pe0) then
4321  active = .true.
4322  call add_pe0_divv(qdt,dtfactor,ixi^l,ixo^l,wct,w,x)
4323  end if
4324  end if
4325 
4326  ! Source for B0 splitting
4327  if (b0field) then
4328  active = .true.
4329  call add_source_b0split(qdt,dtfactor,ixi^l,ixo^l,wct,w,x)
4330  end if
4331 
4332  ! Sources for resistivity in eqs. for e, B1, B2 and B3
4333  if (abs(mhd_eta)>smalldouble)then
4334  active = .true.
4335  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
4336  end if
4337 
4338  if (mhd_eta_hyper>0.d0)then
4339  active = .true.
4340  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
4341  end if
4342 
4343  if(mhd_hydrodynamic_e) then
4344  ! Source for solving hydrodynamic energy
4345  active = .true.
4346  call add_source_hydrodynamic_e(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
4347  else if (mhd_semirelativistic) then
4348  ! add sources for semirelativistic MHD
4349  active = .true.
4350  call add_source_semirelativistic(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
4351  end if
4352  end if
4353 
4354  {^nooned
4355  if(source_split_divb .eqv. qsourcesplit) then
4356  ! Sources related to div B
4357  select case (type_divb)
4358  case (divb_ct)
4359  continue ! Do nothing
4360  case (divb_linde)
4361  active = .true.
4362  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4363  case (divb_glm)
4364  active = .true.
4365  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
4366  case (divb_powel)
4367  active = .true.
4368  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
4369  case (divb_janhunen)
4370  active = .true.
4371  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
4372  case (divb_lindejanhunen)
4373  active = .true.
4374  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4375  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
4376  case (divb_lindepowel)
4377  active = .true.
4378  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4379  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
4380  case (divb_lindeglm)
4381  active = .true.
4382  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
4383  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
4384  case (divb_multigrid)
4385  continue ! Do nothing
4386  case (divb_none)
4387  ! Do nothing
4388  case default
4389  call mpistop('Unknown divB fix')
4390  end select
4391  end if
4392  }
4393 
4394  if(mhd_radiative_cooling) then
4395  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
4396  w,x,qsourcesplit,active, rc_fl)
4397  end if
4398 
4399  if(mhd_viscosity) then
4400  call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
4401  w,x,mhd_energy,qsourcesplit,active)
4402  end if
4403 
4404  if(mhd_gravity) then
4405  call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
4406  w,x,gravity_energy,gravity_rhov,qsourcesplit,active)
4407  end if
4408 
4409  if (mhd_cak_force) then
4410  call cak_add_source(qdt,ixi^l,ixo^l,wct,w,x,mhd_energy,qsourcesplit,active)
4411  end if
4412 
4413  ! update temperature from new pressure, density, and old ionization degree
4414  if(mhd_partial_ionization) then
4415  if(.not.qsourcesplit) then
4416  active = .true.
4417  call mhd_update_temperature(ixi^l,ixo^l,wct,w,x)
4418  end if
4419  end if
4420 
4421  end subroutine mhd_add_source
4422 
4423  subroutine add_pe0_divv(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4425  use mod_geometry
4426 
4427  integer, intent(in) :: ixI^L, ixO^L
4428  double precision, intent(in) :: qdt,dtfactor
4429  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4430  double precision, intent(inout) :: w(ixI^S,1:nw)
4431  double precision :: v(ixI^S,1:ndir)
4432  double precision :: divv(ixI^S)
4433 
4434  call mhd_get_v(wct,x,ixi^l,ixi^l,v)
4435 
4436  if(slab_uniform) then
4437  if(nghostcells .gt. 2) then
4438  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
4439  else
4440  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
4441  end if
4442  else
4443  call divvector(v,ixi^l,ixo^l,divv)
4444  end if
4445  if(local_timestep) then
4446  w(ixo^s,e_)=w(ixo^s,e_)-dtfactor*block%dt(ixo^s)*block%equi_vars(ixo^s,equi_pe0_,0)*divv(ixo^s)
4447  else
4448  w(ixo^s,e_)=w(ixo^s,e_)-qdt*block%equi_vars(ixo^s,equi_pe0_,0)*divv(ixo^s)
4449  endif
4450  end subroutine add_pe0_divv
4451 
4452  !> Compute the Lorentz force (JxB)
4453  subroutine get_lorentz_force(ixI^L,ixO^L,w,JxB)
4455  integer, intent(in) :: ixI^L, ixO^L
4456  double precision, intent(in) :: w(ixI^S,1:nw)
4457  double precision, intent(inout) :: JxB(ixI^S,3)
4458  double precision :: a(ixI^S,3), b(ixI^S,3)
4459  integer :: idir, idirmin
4460  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4461  double precision :: current(ixI^S,7-2*ndir:3)
4462 
4463  b=0.0d0
4464  do idir = 1, ndir
4465  b(ixo^s, idir) = mhd_mag_i_all(w, ixi^l, ixo^l,idir)
4466  end do
4467 
4468  ! store J current in a
4469  call get_current(w,ixi^l,ixo^l,idirmin,current)
4470 
4471  a=0.0d0
4472  do idir=7-2*ndir,3
4473  a(ixo^s,idir)=current(ixo^s,idir)
4474  end do
4475 
4476  call cross_product(ixi^l,ixo^l,a,b,jxb)
4477  end subroutine get_lorentz_force
4478 
4479  !> Compute 1/(1+v_A^2/c^2) for semirelativistic MHD, where v_A is the Alfven
4480  !> velocity
4481  subroutine mhd_gamma2_alfven(ixI^L, ixO^L, w, gamma_A2)
4483  integer, intent(in) :: ixI^L, ixO^L
4484  double precision, intent(in) :: w(ixI^S, nw)
4485  double precision, intent(out) :: gamma_A2(ixO^S)
4486  double precision :: rho(ixI^S)
4487 
4488  ! mhd_get_rho cannot be used as x is not a param
4489  if(has_equi_rho0) then
4490  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
4491  else
4492  rho(ixo^s) = w(ixo^s,rho_)
4493  endif
4494  ! Compute the inverse of 1 + B^2/(rho * c^2)
4495  gamma_a2(ixo^s) = 1.0d0/(1.0d0+mhd_mag_en_all(w, ixi^l, ixo^l)/rho(ixo^s)*inv_squared_c)
4496  end subroutine mhd_gamma2_alfven
4497 
4498  !> Compute 1/sqrt(1+v_A^2/c^2) for semirelativisitic MHD, where v_A is the
4499  !> Alfven velocity
4500  function mhd_gamma_alfven(w, ixI^L, ixO^L) result(gamma_A)
4502  integer, intent(in) :: ixi^l, ixo^l
4503  double precision, intent(in) :: w(ixi^s, nw)
4504  double precision :: gamma_a(ixo^s)
4505 
4506  call mhd_gamma2_alfven(ixi^l, ixo^l, w, gamma_a)
4507  gamma_a = sqrt(gamma_a)
4508  end function mhd_gamma_alfven
4509 
4510  subroutine mhd_get_rho(w,x,ixI^L,ixO^L,rho)
4512  integer, intent(in) :: ixi^l, ixo^l
4513  double precision, intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
4514  double precision, intent(out) :: rho(ixi^s)
4515 
4516  if(has_equi_rho0) then
4517  rho(ixo^s) = w(ixo^s,rho_) + block%equi_vars(ixo^s,equi_rho0_,b0i)
4518  else
4519  rho(ixo^s) = w(ixo^s,rho_)
4520  endif
4521 
4522  end subroutine mhd_get_rho
4523 
4524  !> handle small or negative internal energy
4525  subroutine mhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
4527  use mod_small_values
4528  integer, intent(in) :: ixI^L,ixO^L, ie
4529  double precision, intent(inout) :: w(ixI^S,1:nw)
4530  double precision, intent(in) :: x(ixI^S,1:ndim)
4531  character(len=*), intent(in) :: subname
4532 
4533  integer :: idir
4534  logical :: flag(ixI^S,1:nw)
4535  double precision :: rho(ixI^S)
4536 
4537  flag=.false.
4538  if(has_equi_pe0) then
4539  where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1<small_e)&
4540  flag(ixo^s,ie)=.true.
4541  else
4542  where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
4543  endif
4544  if(any(flag(ixo^s,ie))) then
4545  select case (small_values_method)
4546  case ("replace")
4547  if(has_equi_pe0) then
4548  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
4549  block%equi_vars(ixo^s,equi_pe0_,0)*inv_gamma_1
4550  else
4551  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
4552  endif
4553  case ("average")
4554  call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
4555  case default
4556  ! small values error shows primitive variables
4557  w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
4558  call mhd_get_rho(w,x,ixi^l,ixo^l,rho)
4559  do idir = 1, ndir
4560  w(ixo^s, mom(idir)) = w(ixo^s, mom(idir))/rho(ixo^s)
4561  end do
4562  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
4563  end select
4564  end if
4565 
4566  end subroutine mhd_handle_small_ei
4567 
4568  subroutine mhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
4571 
4572  integer, intent(in) :: ixI^L, ixO^L
4573  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4574  double precision, intent(inout) :: w(ixI^S,1:nw)
4575 
4576  double precision :: iz_H(ixO^S),iz_He(ixO^S), pth(ixI^S)
4577 
4578  call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
4579 
4580  call mhd_get_pthermal(w,x,ixi^l,ixo^l,pth)
4581 
4582  w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
4583  he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
4584 
4585  end subroutine mhd_update_temperature
4586 
4587  !> Source terms after split off time-independent magnetic field
4588  subroutine add_source_b0split(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4590 
4591  integer, intent(in) :: ixI^L, ixO^L
4592  double precision, intent(in) :: qdt, dtfactor,wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4593  double precision, intent(inout) :: w(ixI^S,1:nw)
4594 
4595  double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
4596  integer :: idir
4597 
4598  a=0.d0
4599  b=0.d0
4600  ! for force-free field J0xB0 =0
4601  if(.not.b0field_forcefree) then
4602  ! store B0 magnetic field in b
4603  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
4604 
4605  ! store J0 current in a
4606  do idir=7-2*ndir,3
4607  a(ixo^s,idir)=block%J0(ixo^s,idir)
4608  end do
4609  call cross_product(ixi^l,ixo^l,a,b,axb)
4610  if(local_timestep) then
4611  do idir=1,3
4612  axb(ixo^s,idir)=axb(ixo^s,idir)*block%dt(ixo^s)*dtfactor
4613  enddo
4614  else
4615  axb(ixo^s,:)=axb(ixo^s,:)*qdt
4616  endif
4617  ! add J0xB0 source term in momentum equations
4618  w(ixo^s,mom(1:ndir))=w(ixo^s,mom(1:ndir))+axb(ixo^s,1:ndir)
4619  end if
4620 
4621  if(total_energy) then
4622  a=0.d0
4623  ! for free-free field -(vxB0) dot J0 =0
4624  b(ixo^s,:)=wct(ixo^s,mag(:))
4625  ! store full magnetic field B0+B1 in b
4626  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
4627  ! store velocity in a
4628  call mhd_get_v(wct,x,ixi^l,ixo^l,a(ixi^s,1:ndir))
4629  call cross_product(ixi^l,ixo^l,a,b,axb)
4630  if(local_timestep) then
4631  do idir=1,3
4632  axb(ixo^s,idir)=axb(ixo^s,idir)*block%dt(ixo^s)*dtfactor
4633  enddo
4634  else
4635  axb(ixo^s,:)=axb(ixo^s,:)*qdt
4636  endif
4637  ! add -(vxB) dot J0 source term in energy equation
4638  do idir=7-2*ndir,3
4639  w(ixo^s,e_)=w(ixo^s,e_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
4640  end do
4641  if(mhd_ambipolar) then
4642  !reuse axb
4643  call mhd_get_jxbxb(wct,x,ixi^l,ixo^l,axb)
4644  ! source J0 * E
4645  do idir=sdim,3
4646  !set electric field in jxbxb: E=nuA * jxbxb, where nuA=-etaA/rho^2
4647  call multiplyambicoef(ixi^l,ixo^l,axb(ixi^s,idir),wct,x)
4648  w(ixo^s,e_)=w(ixo^s,e_)+axb(ixo^s,idir)*block%J0(ixo^s,idir)
4649  enddo
4650  endif
4651  end if
4652 
4653  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
4654 
4655  end subroutine add_source_b0split
4656 
4657  !> Source terms for semirelativistic MHD Gombosi 2002 JCP 177, 176
4658  subroutine add_source_semirelativistic(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4660  use mod_geometry
4661 
4662  integer, intent(in) :: ixI^L, ixO^L
4663  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4664  double precision, intent(inout) :: w(ixI^S,1:nw)
4665  double precision, intent(in), optional :: wCTprim(ixI^S,1:nw)
4666 
4667  double precision :: B(ixI^S,3), v(ixI^S,3), E(ixI^S,3), divE(ixI^S)
4668  integer :: idir, idirmin
4669 
4670  ! store B0 magnetic field in b
4671  b=0.d0
4672  if(b0field) then
4673  b(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
4674  else
4675  b(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4676  end if
4677  {^nothreed
4678  v=0.d0
4679  v(ixi^s,1:ndir)=wctprim(ixi^s,mom(1:ndir))
4680 
4681  call cross_product(ixi^l,ixi^l,b,v,e)
4682  }
4683  {^ifthreed
4684  call cross_product(ixi^l,ixi^l,b,wctprim(ixi^s,mom(1:ndir)),e)
4685  }
4686  call divvector(e,ixi^l,ixo^l,dive)
4687  ! curl E => B
4688  call curlvector(e,ixi^l,ixo^l,b,idirmin,1,3)
4689  ! E x (curl E) => v
4690  call cross_product(ixi^l,ixo^l,e,b,v)
4691  ! add source term in momentum equations (1/c0^2-1/c^2)(E dot divE - E x curlE)
4692  ! equation (26) and (27)
4693  do idir=1,ndir
4694  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+qdt*(inv_squared_c0-inv_squared_c)*&
4695  (e(ixo^s,idir)*dive(ixo^s)-v(ixo^s,idir))
4696  end do
4697 
4698  end subroutine add_source_semirelativistic
4699 
4700  !> Source terms for internal energy version of MHD
4701  subroutine add_source_internal_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4703  use mod_geometry
4704 
4705  integer, intent(in) :: ixI^L, ixO^L
4706  double precision, intent(in) :: qdt
4707  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4708  double precision, intent(inout) :: w(ixI^S,1:nw)
4709  double precision, intent(in) :: wCTprim(ixI^S,1:nw)
4710 
4711  double precision :: divv(ixI^S)
4712 
4713  if(slab_uniform) then
4714  if(nghostcells .gt. 2) then
4715  call divvector(wctprim(ixi^s,mom(:)),ixi^l,ixo^l,divv,sixthorder=.true.)
4716  else
4717  call divvector(wctprim(ixi^s,mom(:)),ixi^l,ixo^l,divv,fourthorder=.true.)
4718  end if
4719  else
4720  call divvector(wctprim(ixi^s,mom(:)),ixi^l,ixo^l,divv)
4721  end if
4722  w(ixo^s,e_)=w(ixo^s,e_)-qdt*wctprim(ixo^s,p_)*divv(ixo^s)
4723  if(mhd_ambipolar)then
4724  call add_source_ambipolar_internal_energy(qdt,ixi^l,ixo^l,wct,w,x,e_)
4725  end if
4726 
4727  if(fix_small_values) then
4728  call mhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,'add_source_internal_e')
4729  end if
4730  end subroutine add_source_internal_e
4731 
4732  !> Source terms for hydrodynamic energy version of MHD
4733  subroutine add_source_hydrodynamic_e(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
4735  use mod_geometry
4736  use mod_usr_methods, only: usr_gravity
4737 
4738  integer, intent(in) :: ixI^L, ixO^L
4739  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4740  double precision, intent(inout) :: w(ixI^S,1:nw)
4741  double precision, intent(in), optional :: wCTprim(ixI^S,1:nw)
4742 
4743  double precision :: B(ixI^S,3), J(ixI^S,3), JxB(ixI^S,3)
4744  integer :: idir, idirmin, idims, ix^D
4745  double precision :: current(ixI^S,7-2*ndir:3)
4746  double precision :: bu(ixO^S,1:ndir), tmp(ixO^S), b2(ixO^S)
4747  double precision :: gravity_field(ixI^S,1:ndir), Vaoc
4748 
4749  {^nothreed
4750  b=0.0d0
4751  do idir = 1, ndir
4752  b(ixo^s, idir) = wct(ixo^s,mag(idir))
4753  end do
4754 
4755  !call get_current(wCT,ixI^L,ixO^L,idirmin,current)
4756  call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,7-2*ndir,ndir,.true.)
4757 
4758  j=0.0d0
4759  do idir=7-2*ndir,3
4760  j(ixo^s,idir)=current(ixo^s,idir)
4761  end do
4762 
4763  ! get Lorentz force JxB
4764  call cross_product(ixi^l,ixo^l,j,b,jxb)
4765  }
4766  {^ifthreed
4767  !call get_current(wCT,ixI^L,ixO^L,idirmin,current)
4768  ! get current in fourth order accuracy in Cartesian
4769  call curlvector(wct(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,1,ndir,.true.)
4770  ! get Lorentz force JxB
4771  call cross_product(ixi^l,ixo^l,current,wct(ixi^s,mag(1:ndir)),jxb)
4772  }
4773 
4774  if(mhd_semirelativistic) then
4775  ! (v . nabla) v
4776  do idir=1,ndir
4777  do idims=1,ndim
4778  call gradient(wctprim(ixi^s,mom(idir)),ixi^l,ixo^l,idims,j(ixi^s,idims))
4779  end do
4780  b(ixo^s,idir)=sum(wctprim(ixo^s,mom(1:ndir))*j(ixo^s,1:ndir),dim=ndim+1)
4781  end do
4782  ! nabla p
4783  do idir=1,ndir
4784  call gradient(wctprim(ixi^s,p_),ixi^l,ixo^l,idir,j(ixi^s,idir))
4785  end do
4786 
4787  if(mhd_gravity) then
4788  gravity_field=0.d0
4789  call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field(ixi^s,1:ndim))
4790  do idir=1,ndir
4791  b(ixo^s,idir)=wct(ixo^s,rho_)*(b(ixo^s,idir)-gravity_field(ixo^s,idir))+j(ixo^s,idir)-jxb(ixo^s,idir)
4792  end do
4793  else
4794  do idir=1,ndir
4795  b(ixo^s,idir)=wct(ixo^s,rho_)*b(ixo^s,idir)+j(ixo^s,idir)-jxb(ixo^s,idir)
4796  end do
4797  end if
4798 
4799  b2(ixo^s)=sum(wct(ixo^s,mag(:))**2,dim=ndim+1)
4800  tmp(ixo^s)=sqrt(b2(ixo^s))
4801  where(tmp(ixo^s)>smalldouble)
4802  tmp(ixo^s)=1.d0/tmp(ixo^s)
4803  else where
4804  tmp(ixo^s)=0.d0
4805  end where
4806  ! unit vector of magnetic field
4807  do idir=1,ndir
4808  bu(ixo^s,idir)=wct(ixo^s,mag(idir))*tmp(ixo^s)
4809  end do
4810 
4811  !b2(ixO^S)=b2(ixO^S)/w(ixO^S,rho_)*inv_squared_c
4812  !b2(ixO^S)=b2(ixO^S)/(1.d0+b2(ixO^S))
4813  {do ix^db=ixomin^db,ixomax^db\}
4814  ! Va^2/c^2
4815  vaoc=b2(ix^d)/w(ix^d,rho_)*inv_squared_c
4816  ! Va^2/c^2 / (1+Va^2/c^2)
4817  b2(ix^d)=vaoc/(1.d0+vaoc)
4818  {end do\}
4819  ! bu . F
4820  tmp(ixo^s)=sum(bu(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
4821  ! Rempel 2017 ApJ 834, 10 equation (54)
4822  do idir=1,ndir
4823  j(ixo^s,idir)=b2(ixo^s)*(b(ixo^s,idir)-bu(ixo^s,idir)*tmp(ixo^s))
4824  end do
4825  !! Rempel 2017 ApJ 834, 10 equation (29) add SR force at momentum equation
4826  do idir=1,ndir
4827  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+qdt*j(ixo^s,idir)
4828  end do
4829  ! Rempel 2017 ApJ 834, 10 equation (30) add work of Lorentz force and SR force
4830  w(ixo^s,e_)=w(ixo^s,e_)+qdt*sum(wctprim(ixo^s,mom(1:ndir))*&
4831  (jxb(ixo^s,1:ndir)+j(ixo^s,1:ndir)),dim=ndim+1)
4832  else
4833  ! add work of Lorentz force
4834  w(ixo^s,e_)=w(ixo^s,e_)+qdt*sum(wctprim(ixo^s,mom(1:ndir))*jxb(ixo^s,1:ndir),dim=ndim+1)
4835  end if
4836 
4837  end subroutine add_source_hydrodynamic_e
4838 
4839  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
4840  !> each direction, non-conservative. If the fourthorder precompiler flag is
4841  !> set, uses fourth order central difference for the laplacian. Then the
4842  !> stencil is 5 (2 neighbours).
4843  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
4845  use mod_usr_methods
4846  use mod_geometry
4847 
4848  integer, intent(in) :: ixI^L, ixO^L
4849  double precision, intent(in) :: qdt
4850  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4851  double precision, intent(inout) :: w(ixI^S,1:nw)
4852  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
4853  integer :: lxO^L, kxO^L
4854 
4855  double precision :: tmp(ixI^S),tmp2(ixI^S)
4856 
4857  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4858  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4859  double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
4860 
4861  ! Calculating resistive sources involve one extra layer
4862  if (mhd_4th_order) then
4863  ixa^l=ixo^l^ladd2;
4864  else
4865  ixa^l=ixo^l^ladd1;
4866  end if
4867 
4868  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4869  call mpistop("Error in add_source_res1: Non-conforming input limits")
4870 
4871  ! Calculate current density and idirmin
4872  call get_current(wct,ixi^l,ixo^l,idirmin,current)
4873 
4874  if (mhd_eta>zero)then
4875  eta(ixa^s)=mhd_eta
4876  gradeta(ixo^s,1:ndim)=zero
4877  else
4878  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
4879  ! assumes that eta is not function of current?
4880  do idim=1,ndim
4881  call gradient(eta,ixi^l,ixo^l,idim,tmp)
4882  gradeta(ixo^s,idim)=tmp(ixo^s)
4883  end do
4884  end if
4885 
4886  if(b0field) then
4887  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
4888  else
4889  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
4890  end if
4891 
4892  do idir=1,ndir
4893  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
4894  if (mhd_4th_order) then
4895  tmp(ixo^s)=zero
4896  tmp2(ixi^s)=bf(ixi^s,idir)
4897  do idim=1,ndim
4898  lxo^l=ixo^l+2*kr(idim,^d);
4899  jxo^l=ixo^l+kr(idim,^d);
4900  hxo^l=ixo^l-kr(idim,^d);
4901  kxo^l=ixo^l-2*kr(idim,^d);
4902  tmp(ixo^s)=tmp(ixo^s)+&
4903  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
4904  /(12.0d0 * dxlevel(idim)**2)
4905  end do
4906  else
4907  tmp(ixo^s)=zero
4908  tmp2(ixi^s)=bf(ixi^s,idir)
4909  do idim=1,ndim
4910  jxo^l=ixo^l+kr(idim,^d);
4911  hxo^l=ixo^l-kr(idim,^d);
4912  tmp(ixo^s)=tmp(ixo^s)+&
4913  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
4914  end do
4915  end if
4916 
4917  ! Multiply by eta
4918  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
4919 
4920  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
4921  if (mhd_eta<zero)then
4922  do jdir=1,ndim; do kdir=idirmin,3
4923  if (lvc(idir,jdir,kdir)/=0)then
4924  if (lvc(idir,jdir,kdir)==1)then
4925  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4926  else
4927  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
4928  end if
4929  end if
4930  end do; end do
4931  end if
4932 
4933  ! Add sources related to eta*laplB-grad(eta) x J to B and e
4934  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
4935  if(total_energy) then
4936  w(ixo^s,e_)=w(ixo^s,e_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
4937  end if
4938  end do ! idir
4939 
4940  if(mhd_energy) then
4941  ! de/dt+=eta*J**2
4942  w(ixo^s,e_)=w(ixo^s,e_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4943  end if
4944 
4945  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
4946 
4947  end subroutine add_source_res1
4948 
4949  !> Add resistive source to w within ixO
4950  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
4951  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4953  use mod_usr_methods
4954  use mod_geometry
4955 
4956  integer, intent(in) :: ixI^L, ixO^L
4957  double precision, intent(in) :: qdt
4958  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4959  double precision, intent(inout) :: w(ixI^S,1:nw)
4960 
4961  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4962  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4963  double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4964  integer :: ixA^L,idir,idirmin,idirmin1
4965 
4966  ixa^l=ixo^l^ladd2;
4967 
4968  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4969  call mpistop("Error in add_source_res2: Non-conforming input limits")
4970 
4971  ixa^l=ixo^l^ladd1;
4972  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
4973  ! Determine exact value of idirmin while doing the loop.
4974  call get_current(wct,ixi^l,ixa^l,idirmin,current)
4975 
4976  tmpvec=zero
4977  if(mhd_eta>zero)then
4978  do idir=idirmin,3
4979  tmpvec(ixa^s,idir)=current(ixa^s,idir)*mhd_eta
4980  end do
4981  else
4982  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
4983  do idir=idirmin,3
4984  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4985  end do
4986  end if
4987 
4988  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
4989  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4990  if(stagger_grid) then
4991  if(ndim==2.and.ndir==3) then
4992  ! if 2.5D
4993  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4994  end if
4995  else
4996  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4997  end if
4998 
4999  if(mhd_energy) then
5000  if(mhd_eta>zero)then
5001  tmp(ixo^s)=qdt*mhd_eta*sum(current(ixo^s,:)**2,dim=ndim+1)
5002  else
5003  tmp(ixo^s)=qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
5004  end if
5005  if(total_energy) then
5006  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
5007  ! de1/dt= eta J^2 - B1 dot curl(eta J)
5008  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)-&
5009  qdt*sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1)
5010  else
5011  ! add eta*J**2 source term in the internal energy equation
5012  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)
5013  end if
5014  end if
5015 
5016  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
5017  end subroutine add_source_res2
5018 
5019  !> Add Hyper-resistive source to w within ixO
5020  !> Uses 9 point stencil (4 neighbours) in each direction.
5021  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
5023  use mod_geometry
5024 
5025  integer, intent(in) :: ixI^L, ixO^L
5026  double precision, intent(in) :: qdt
5027  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5028  double precision, intent(inout) :: w(ixI^S,1:nw)
5029  !.. local ..
5030  double precision :: current(ixI^S,7-2*ndir:3)
5031  double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
5032  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
5033 
5034  ixa^l=ixo^l^ladd3;
5035  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
5036  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
5037 
5038  call get_current(wct,ixi^l,ixa^l,idirmin,current)
5039  tmpvec(ixa^s,1:ndir)=zero
5040  do jdir=idirmin,3
5041  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
5042  end do
5043 
5044  ixa^l=ixo^l^ladd2;
5045  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5046 
5047  ixa^l=ixo^l^ladd1;
5048  tmpvec(ixa^s,1:ndir)=zero
5049  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
5050  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mhd_eta_hyper
5051 
5052  ixa^l=ixo^l;
5053  tmpvec2(ixa^s,1:ndir)=zero
5054  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
5055 
5056  do idir=1,ndir
5057  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
5058  end do
5059 
5060  if(total_energy) then
5061  ! de/dt= +div(B x Ehyper)
5062  ixa^l=ixo^l^ladd1;
5063  tmpvec2(ixa^s,1:ndir)=zero
5064  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
5065  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
5066  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
5067  end do; end do; end do
5068  tmp(ixo^s)=zero
5069  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
5070  w(ixo^s,e_)=w(ixo^s,e_)+tmp(ixo^s)*qdt
5071  end if
5072 
5073  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
5074 
5075  end subroutine add_source_hyperres
5076 
5077  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
5078  ! Add divB related sources to w within ixO
5079  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
5080  ! giving the EGLM-MHD scheme or GLM-MHD scheme
5082  use mod_geometry
5083 
5084  integer, intent(in) :: ixI^L, ixO^L
5085  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5086  double precision, intent(inout) :: w(ixI^S,1:nw)
5087  double precision:: divb(ixI^S)
5088  integer :: idim,idir
5089  double precision :: gradPsi(ixI^S)
5090 
5091 
5092  ! dPsi/dt = - Ch^2/Cp^2 Psi
5093  if (mhd_glm_alpha < zero) then
5094  w(ixo^s,psi_) = abs(mhd_glm_alpha)*wct(ixo^s,psi_)
5095  else
5096  ! implicit update of Psi variable
5097  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
5098  if(slab_uniform) then
5099  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
5100  else
5101  w(ixo^s,psi_) = dexp(-qdt*cmax_global*mhd_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
5102  end if
5103  end if
5104 
5105  if(mhd_glm_extended) then
5106  ! gradient of Psi
5107  if(total_energy) then
5108  do idim=1,ndim
5109  select case(typegrad)
5110  case("central")
5111  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
5112  case("limited")
5113  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
5114  end select
5115  ! e = e -qdt (b . grad(Psi))
5116  w(ixo^s,e_) = w(ixo^s,e_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
5117  end do
5118  end if
5119 
5120  ! We calculate now div B
5121  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
5122 
5123  ! m = m - qdt b div b
5124  do idir=1,ndir
5125  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
5126  end do
5127  end if
5128 
5129  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
5130 
5131  end subroutine add_source_glm
5132 
5133  !> Add divB related sources to w within ixO corresponding to Powel
5134  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
5136 
5137  integer, intent(in) :: ixI^L, ixO^L
5138  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5139  double precision, intent(inout) :: w(ixI^S,1:nw)
5140  double precision :: divb(ixI^S),v(ixI^S,1:ndir)
5141  integer :: idir
5142 
5143  ! We calculate now div B
5144  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
5145 
5146  ! calculate velocity
5147  call mhd_get_v(wct,x,ixi^l,ixo^l,v)
5148 
5149  if (total_energy) then
5150  ! e = e - qdt (v . b) * div b
5151  w(ixo^s,e_)=w(ixo^s,e_)-&
5152  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
5153  end if
5154 
5155  ! b = b - qdt v * div b
5156  do idir=1,ndir
5157  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
5158  end do
5159 
5160  ! m = m - qdt b div b
5161  do idir=1,ndir
5162  w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-qdt*mhd_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
5163  end do
5164 
5165  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
5166 
5167  end subroutine add_source_powel
5168 
5169  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
5170  ! Add divB related sources to w within ixO
5171  ! corresponding to Janhunen, just the term in the induction equation.
5173 
5174  integer, intent(in) :: ixI^L, ixO^L
5175  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5176  double precision, intent(inout) :: w(ixI^S,1:nw)
5177  double precision :: divb(ixI^S),vel(ixI^S,1:ndir)
5178  integer :: idir
5179 
5180  ! We calculate now div B
5181  call get_divb(wct,ixi^l,ixo^l,divb, mhd_divb_4thorder)
5182 
5183  call mhd_get_v(wct,x,ixi^l,ixo^l,vel)
5184  ! b = b - qdt v * div b
5185  do idir=1,ndir
5186  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s,idir)*divb(ixo^s)
5187  end do
5188 
5189  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
5190 
5191  end subroutine add_source_janhunen
5192 
5193  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
5194  ! Add Linde's divB related sources to wnew within ixO
5196  use mod_geometry
5197 
5198  integer, intent(in) :: ixI^L, ixO^L
5199  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
5200  double precision, intent(inout) :: w(ixI^S,1:nw)
5201  integer :: idim, idir, ixp^L, i^D, iside
5202  double precision :: divb(ixI^S),graddivb(ixI^S)
5203  logical, dimension(-1:1^D&) :: leveljump
5204 
5205  ! Calculate div B
5206  ixp^l=ixo^l^ladd1;
5207  call get_divb(wct,ixi^l,ixp^l,divb, mhd_divb_4thorder)
5208 
5209  ! for AMR stability, retreat one cell layer from the boarders of level jump
5210  {do i^db=-1,1\}
5211  if(i^d==0|.and.) cycle
5212  if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
5213  leveljump(i^d)=.true.
5214  else
5215  leveljump(i^d)=.false.
5216  end if
5217  {end do\}
5218 
5219  ixp^l=ixo^l;
5220  do idim=1,ndim
5221  select case(idim)
5222  {case(^d)
5223  do iside=1,2
5224  i^dd=kr(^dd,^d)*(2*iside-3);
5225  if (leveljump(i^dd)) then
5226  if (iside==1) then
5227  ixpmin^d=ixomin^d-i^d
5228  else
5229  ixpmax^d=ixomax^d-i^d
5230  end if
5231  end if
5232  end do
5233  \}
5234  end select
5235  end do
5236 
5237  ! Add Linde's diffusive terms
5238  do idim=1,ndim
5239  ! Calculate grad_idim(divb)
5240  select case(typegrad)
5241  case("central")
5242  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
5243  case("limited")
5244  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
5245  end select
5246 
5247  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
5248  if (slab_uniform) then
5249  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
5250  else
5251  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
5252  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
5253  end if
5254 
5255  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
5256 
5257  if (typedivbdiff=='all' .and. total_energy) then
5258  ! e += B_idim*eta*grad_idim(divb)
5259  w(ixp^s,e_)=w(ixp^s,e_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
5260  end if
5261  end do
5262 
5263  if (fix_small_values) call mhd_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
5264 
5265  end subroutine add_source_linde
5266 
5267 
5268  !> get dimensionless div B = |divB| * volume / area / |B|
5269  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
5270 
5272 
5273  integer, intent(in) :: ixi^l, ixo^l
5274  double precision, intent(in) :: w(ixi^s,1:nw)
5275  double precision :: divb(ixi^s), dsurface(ixi^s)
5276 
5277  double precision :: invb(ixo^s)
5278  integer :: ixa^l,idims
5279 
5280  call get_divb(w,ixi^l,ixo^l,divb)
5281  invb(ixo^s)=sqrt(mhd_mag_en_all(w,ixi^l,ixo^l))
5282  where(invb(ixo^s)/=0.d0)
5283  invb(ixo^s)=1.d0/invb(ixo^s)
5284  end where
5285  if(slab_uniform) then
5286  divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
5287  else
5288  ixamin^d=ixomin^d-1;
5289  ixamax^d=ixomax^d-1;
5290  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
5291  do idims=1,ndim
5292  ixa^l=ixo^l-kr(idims,^d);
5293  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
5294  end do
5295  divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
5296  block%dvolume(ixo^s)/dsurface(ixo^s)
5297  end if
5298 
5299  end subroutine get_normalized_divb
5300 
5301  !> Calculate idirmin and the idirmin:3 components of the common current array
5302  !> make sure that dxlevel(^D) is set correctly.
5303  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
5305  use mod_geometry
5306 
5307  integer, intent(in) :: ixo^l, ixi^l
5308  double precision, intent(in) :: w(ixi^s,1:nw)
5309  integer, intent(out) :: idirmin
5310  integer :: idir, idirmin0
5311 
5312  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
5313  double precision :: current(ixi^s,7-2*ndir:3)
5314 
5315  idirmin0 = 7-2*ndir
5316 
5317  call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
5318 
5319  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
5320  block%J0(ixo^s,idirmin0:3)
5321  end subroutine get_current
5322 
5323  !> If resistivity is not zero, check diffusion time limit for dt
5324  subroutine mhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
5326  use mod_usr_methods
5328  use mod_viscosity, only: viscosity_get_dt
5329  use mod_gravity, only: gravity_get_dt
5330  use mod_cak_force, only: cak_get_dt
5331 
5332  integer, intent(in) :: ixI^L, ixO^L
5333  double precision, intent(inout) :: dtnew
5334  double precision, intent(in) :: dx^D
5335  double precision, intent(in) :: w(ix