MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_twofl_phys.t
Go to the documentation of this file.
1 !> Magneto-hydrodynamics module
3 
4 #include "amrvac.h"
5 
6  use mod_physics
7  use mod_global_parameters, only: std_len
12  use mod_comm_lib, only: mpistop
13 
14  implicit none
15  private
16  !! E_c = E_kin + E_mag + E_int
17  !! E_n = E_kin + E_int
18  integer, public, parameter :: eq_energy_tot=2
19  !! E_c = E_int
20  !! E_n = E_int
21  integer, public, parameter :: eq_energy_int=1
22  !! E_n, E_c are calculated from density as c_adiab rho^gamma
23  !! No energy equation => no variable assigned for it
24  integer, public, parameter :: eq_energy_none=0
25  !! E_c = E_kin + E_int
26  !! E_n = E_kin + E_int
27  integer, public, parameter :: eq_energy_ki=3
28 
29  integer, public, protected :: twofl_eq_energy = eq_energy_tot
30 
31  !> Whether hyperdiffusivity is used
32  logical, public, protected :: twofl_hyperdiffusivity = .false.
33  logical, public, protected :: twofl_dump_hyperdiffusivity_coef = .false.
34  double precision, public, protected, allocatable :: c_shk(:)
35  double precision, public, protected, allocatable :: c_hyp(:)
36 
37  !> Whether thermal conduction is used
38  logical, public, protected :: twofl_thermal_conduction_c = .false.
39  !> type of TC used: 1: adapted module (mhd implementation), 2: adapted module (hd implementation)
40  integer, parameter, private :: mhd_tc =1
41  integer, parameter, private :: hd_tc =2
42  integer, protected :: use_twofl_tc_c = mhd_tc
43 
44  !> Whether radiative cooling is added
45  logical, public, protected :: twofl_radiative_cooling_c = .false.
46  type(rc_fluid), public, allocatable :: rc_fl_c
47 
48  !> Whether viscosity is added
49  logical, public, protected :: twofl_viscosity = .false.
50 
51  !> Whether gravity is added: common flag for charges and neutrals
52  logical, public, protected :: twofl_gravity = .false.
53 
54  !> whether dump full variables (when splitting is used) in a separate dat file
55  logical, public, protected :: twofl_dump_full_vars = .false.
56 
57  !> Whether Hall-MHD is used
58  logical, public, protected :: twofl_hall = .false.
59 
60  type(tc_fluid), public, allocatable :: tc_fl_c
61  type(te_fluid), public, allocatable :: te_fl_c
62 
63  type(tc_fluid), allocatable :: tc_fl_n
64  logical, public, protected :: twofl_thermal_conduction_n = .false.
65  logical, public, protected :: twofl_radiative_cooling_n = .false.
66  type(rc_fluid), allocatable :: rc_fl_n
67 
68  !> Whether TRAC method is used
69  logical, public, protected :: twofl_trac = .false.
70 
71  !> Whether GLM-MHD is used
72  logical, public, protected :: twofl_glm = .false.
73 
74  !> Which TRAC method is used
75  integer, public, protected :: twofl_trac_type=1
76 
77  !> Height of the mask used in the TRAC method
78  double precision, public, protected :: twofl_trac_mask = 0.d0
79 
80  !> Whether divB cleaning sources are added splitting from fluid solver
81  logical, public, protected :: source_split_divb = .false.
82 
83  !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
84  !> taking values within [0, 1]
85  double precision, public :: twofl_glm_alpha = 0.5d0
86 
87  !> MHD fourth order
88  logical, public, protected :: twofl_4th_order = .false.
89 
90  !> Index of the density (in the w array)
91  integer, public :: rho_c_
92 
93  !> Indices of the momentum density
94  integer, allocatable, public :: mom_c(:)
95 
96  !> Index of the energy density (-1 if not present)
97  integer, public :: e_c_=-1
98 
99  !> Index of the cutoff temperature for the TRAC method
100  integer, public :: tcoff_c_
101  integer, public :: tweight_c_
102 
103  !> Indices of the GLM psi
104  integer, public, protected :: psi_
105 
106  !> equi vars flags
107  logical, public :: has_equi_rho_c0 = .false.
108  logical, public :: has_equi_pe_c0 = .false.
109 
110  !> equi vars indices in the state%equi_vars array
111  integer, public :: equi_rho_c0_ = -1
112  integer, public :: equi_pe_c0_ = -1
113  logical, public :: twofl_equi_thermal_c = .false.
114 
115  logical, public :: twofl_equi_thermal = .false.
116  !neutrals:
117 
118  integer, public :: rho_n_
119  integer, allocatable, public :: mom_n(:)
120  integer, public :: e_n_
121  integer, public :: tcoff_n_
122  integer, public :: tweight_n_
123  logical, public :: has_equi_rho_n0 = .false.
124  logical, public :: has_equi_pe_n0 = .false.
125  integer, public :: equi_rho_n0_ = -1
126  integer, public :: equi_pe_n0_ = -1
127 
128  ! related to collisions:
129  !> collisional alpha
130  double precision, public :: twofl_alpha_coll = 0d0
131  logical, public :: twofl_alpha_coll_constant = .true.
132  !> whether include thermal exchange collisional terms
133  logical, public :: twofl_coll_inc_te = .true.
134  !> whether include ionization/recombination inelastic collisional terms
135  logical, public :: twofl_coll_inc_ionrec = .false.
136  logical, public :: twofl_equi_thermal_n = .false.
137  double precision, public :: dtcollpar = -1d0 !negative value does not impose restriction on the timestep
138  !> whether dump collisional terms in a separte dat file
139  logical, public, protected :: twofl_dump_coll_terms = .false.
140 
141  ! TODO Helium abundance not used, radiative cooling init uses it
142  ! not in parameters list anymore
143  double precision, public, protected :: he_abundance = 0d0
144  ! two fluid is only H plasma
145  double precision, public, protected :: rc = 2d0
146  double precision, public, protected :: rn = 1d0
147 
148  !> The adiabatic index
149  double precision, public :: twofl_gamma = 5.d0/3.0d0
150 
151  !> The adiabatic constant
152  double precision, public :: twofl_adiab = 1.0d0
153 
154  !> The MHD resistivity
155  double precision, public :: twofl_eta = 0.0d0
156 
157  !> The MHD hyper-resistivity
158  double precision, public :: twofl_eta_hyper = 0.0d0
159 
160  !> The MHD Hall coefficient
161  double precision, public :: twofl_etah = 0.0d0
162 
163  !> The small_est allowed energy
164  double precision, protected :: small_e
165 
166  !> Method type to clean divergence of B
167  character(len=std_len), public, protected :: typedivbfix = 'linde'
168 
169  !> Method type of constrained transport
170  character(len=std_len), public, protected :: type_ct = 'uct_contact'
171 
172  !> Whether divB is computed with a fourth order approximation
173  logical, public, protected :: twofl_divb_4thorder = .false.
174 
175  !> Method type in a integer for good performance
176  integer :: type_divb
177 
178  !> Coefficient of diffusive divB cleaning
179  double precision :: divbdiff = 0.8d0
180 
181  !> Update all equations due to divB cleaning
182  character(len=std_len) :: typedivbdiff = 'all'
183 
184  !> clean initial divB
185  logical, public :: clean_initial_divb = .false.
186 
187  !> Add divB wave in Roe solver
188  logical, public :: divbwave = .true.
189 
190  !> To control divB=0 fix for boundary
191  logical, public, protected :: boundary_divbfix(2*^nd)=.true.
192 
193  !> To skip * layer of ghost cells during divB=0 fix for boundary
194  integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
195 
196  !> B0 field is force-free
197  logical, public, protected :: b0field_forcefree=.true.
198 
199  logical :: twofl_cbounds_species = .true.
200 
201  !> added from modules: gravity
202  !> source split or not
203  logical :: grav_split= .false.
204 
205  !> gamma minus one and its inverse
206  double precision :: gamma_1, inv_gamma_1
207 
208  ! DivB cleaning methods
209  integer, parameter :: divb_none = 0
210  integer, parameter :: divb_multigrid = -1
211  integer, parameter :: divb_glm = 1
212  integer, parameter :: divb_powel = 2
213  integer, parameter :: divb_janhunen = 3
214  integer, parameter :: divb_linde = 4
215  integer, parameter :: divb_lindejanhunen = 5
216  integer, parameter :: divb_lindepowel = 6
217  integer, parameter :: divb_lindeglm = 7
218  integer, parameter :: divb_ct = 8
219 
220  ! Public methods
221  public :: twofl_phys_init
222  public :: twofl_to_conserved
223  public :: twofl_to_primitive
224  public :: get_divb
225  public :: get_rhoc_tot
226  public :: twofl_get_v_c_idim
227  ! TODO needed for the roe, see if can be used for n
229  public :: get_rhon_tot
230  public :: get_alpha_coll
231  public :: get_gamma_ion_rec
232  public :: twofl_get_v_n_idim
233  public :: get_current
234  public :: twofl_get_pthermal_c
235  public :: twofl_get_pthermal_n
236  public :: twofl_face_to_center
237  public :: get_normalized_divb
238  public :: b_from_vector_potential
239  public :: usr_mask_gamma_ion_rec
240  public :: usr_mask_alpha
241 
242  {^nooned
244  }
245 
246  abstract interface
247 
248  subroutine implicit_mult_factor_subroutine(ixI^L, ixO^L, step_dt, JJ, res)
249  integer, intent(in) :: ixi^l, ixo^l
250  double precision, intent(in) :: step_dt
251  double precision, intent(in) :: jj(ixi^s)
252  double precision, intent(out) :: res(ixi^s)
253 
254  end subroutine implicit_mult_factor_subroutine
255 
256  subroutine mask_subroutine(ixI^L,ixO^L,w,x,res)
258  integer, intent(in) :: ixi^l, ixo^l
259  double precision, intent(in) :: x(ixi^s,1:ndim)
260  double precision, intent(in) :: w(ixi^s,1:nw)
261  double precision, intent(inout) :: res(ixi^s)
262  end subroutine mask_subroutine
263 
264  subroutine mask_subroutine2(ixI^L,ixO^L,w,x,res1, res2)
266  integer, intent(in) :: ixI^L, ixO^L
267  double precision, intent(in) :: x(ixI^S,1:ndim)
268  double precision, intent(in) :: w(ixI^S,1:nw)
269  double precision, intent(inout) :: res1(ixI^S),res2(ixI^S)
270  end subroutine mask_subroutine2
271 
272  end interface
273 
274  procedure(implicit_mult_factor_subroutine), pointer :: calc_mult_factor => null()
275  integer, protected :: twofl_implicit_calc_mult_method = 1
276  procedure(mask_subroutine), pointer :: usr_mask_alpha => null()
277  procedure(mask_subroutine2), pointer :: usr_mask_gamma_ion_rec => null()
278 
279 contains
280 
281  !> Read this module"s parameters from a file
282  subroutine twofl_read_params(files)
284  character(len=*), intent(in) :: files(:)
285  integer :: n
286 
287  namelist /twofl_list/ twofl_eq_energy, twofl_gamma, twofl_adiab,&
291  typedivbdiff, type_ct, divbwave, si_unit, b0field,&
298  twofl_dump_coll_terms,twofl_implicit_calc_mult_method,&
301  twofl_trac, twofl_trac_type, twofl_trac_mask,twofl_cbounds_species
302 
303  do n = 1, size(files)
304  open(unitpar, file=trim(files(n)), status="old")
305  read(unitpar, twofl_list, end=111)
306 111 close(unitpar)
307  end do
308 
309  end subroutine twofl_read_params
310 
311  subroutine twofl_init_hyper(files)
314  character(len=*), intent(in) :: files(:)
315  integer :: n
316 
317  namelist /hyperdiffusivity_list/ c_shk, c_hyp
318 
319  do n = 1, size(files)
320  open(unitpar, file=trim(files(n)), status="old")
321  read(unitpar, hyperdiffusivity_list, end=113)
322 113 close(unitpar)
323  end do
324 
325  call hyperdiffusivity_init()
326 
327  !!DEBUG
328  if(mype .eq. 0) then
329  print*, "Using Hyperdiffusivity"
330  print*, "C_SHK ", c_shk(:)
331  print*, "C_HYP ", c_hyp(:)
332  endif
333 
334  end subroutine twofl_init_hyper
335 
336  !> Write this module's parameters to a snapsoht
337  subroutine twofl_write_info(fh)
339  integer, intent(in) :: fh
340  integer, parameter :: n_par = 1
341  double precision :: values(n_par)
342  character(len=name_len) :: names(n_par)
343  integer, dimension(MPI_STATUS_SIZE) :: st
344  integer :: er
345 
346  call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
347 
348  names(1) = "gamma"
349  values(1) = twofl_gamma
350  call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
351  call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
352  end subroutine twofl_write_info
353 
354  subroutine twofl_phys_init()
358  use mod_viscosity, only: viscosity_init
359  !use mod_gravity, only: gravity_init
362  {^nooned
364  }
365  integer :: itr, idir
366 
367  call twofl_read_params(par_files)
368  physics_type = "twofl"
369  if (twofl_cbounds_species) then
370  number_species = 2
371  endif
372  phys_energy=.true.
373  !> Solve total energy equation or not
374  ! for the two fluid the true value means
375  ! E_charges = E_mag + E_kin_charges + E_int_charges
376  ! E_neutrals = E_kin_neutrals + E_int_neutrals
377  phys_total_energy=.false.
378 
379  !> Solve internal energy instead of total energy
380  ! for the two fluid the true value means
381  ! E_charges = E_int_charges
382  ! E_neutrals = E_int_neutrals
383  phys_internal_e=.false.
384 
385  ! For the two fluid phys_energy=.true. and phys_internal_e=.false. and phys_total_energy = .false. means
386  ! E_charges = E_kin_charges + E_int_charges
387  ! E_neutrals = E_kin_neutrals + E_int_neutrals
388  phys_gamma = twofl_gamma
389 
390  if(twofl_eq_energy == eq_energy_int) then
391  phys_internal_e = .true.
392  elseif(twofl_eq_energy == eq_energy_tot) then
393  phys_total_energy = .true.
394  elseif(twofl_eq_energy == eq_energy_none) then
395  phys_energy = .false.
396  endif
397 
400 
401  if(.not. phys_energy) then
404  if(mype==0) write(*,*) 'WARNING: set twofl_thermal_conduction_n=F when twofl_energy=F'
405  end if
408  if(mype==0) write(*,*) 'WARNING: set twofl_radiative_cooling_n=F when twofl_energy=F'
409  end if
412  if(mype==0) write(*,*) 'WARNING: set twofl_thermal_conduction_c=F when twofl_energy=F'
413  end if
416  if(mype==0) write(*,*) 'WARNING: set twofl_radiative_cooling_c=F when twofl_energy=F'
417  end if
418  if(twofl_trac) then
419  twofl_trac=.false.
420  if(mype==0) write(*,*) 'WARNING: set twofl_trac=F when twofl_energy=F'
421  end if
422  end if
423  {^ifoned
424  if(twofl_trac .and. twofl_trac_type .gt. 1) then
426  if(mype==0) write(*,*) 'WARNING: set twofl_trac_type=1 for 1D simulation'
427  end if
428  }
429  if(twofl_trac .and. twofl_trac_type .le. 3) then
430  twofl_trac_mask=bigdouble
431  if(mype==0) write(*,*) 'WARNING: set twofl_trac_mask==bigdouble for global TRAC method'
432  end if
434 
435  ! set default gamma for polytropic/isothermal process
436  if(ndim==1) typedivbfix='none'
437  select case (typedivbfix)
438  case ('none')
439  type_divb = divb_none
440  {^nooned
441  case ('multigrid')
442  type_divb = divb_multigrid
443  use_multigrid = .true.
444  mg%operator_type = mg_laplacian
445  phys_global_source_after => twofl_clean_divb_multigrid
446  }
447  case ('glm')
448  twofl_glm = .true.
449  need_global_cmax = .true.
450  type_divb = divb_glm
451  case ('powel', 'powell')
452  type_divb = divb_powel
453  case ('janhunen')
454  type_divb = divb_janhunen
455  case ('linde')
456  type_divb = divb_linde
457  case ('lindejanhunen')
458  type_divb = divb_lindejanhunen
459  case ('lindepowel')
460  type_divb = divb_lindepowel
461  case ('lindeglm')
462  twofl_glm = .true.
463  need_global_cmax = .true.
464  type_divb = divb_lindeglm
465  case ('ct')
466  type_divb = divb_ct
467  stagger_grid = .true.
468  case default
469  call mpistop('Unknown divB fix')
470  end select
471 
472  allocate(start_indices(number_species))
473  allocate(stop_indices(number_species))
474  start_indices(1)=1
475  !allocate charges first and the same order as in mhd module
476  rho_c_ = var_set_fluxvar("rho_c", "rho_c")
477  !set variables from mod_variables to point to charges vars
478  iw_rho = rho_c_
479 
480  allocate(mom_c(ndir))
481  do idir=1,ndir
482  mom_c(idir) = var_set_fluxvar("m_c","v_c",idir)
483  enddo
484 
485  allocate(iw_mom(ndir))
486  iw_mom(1:ndir) = mom_c(1:ndir)
487 
488  ! Set index of energy variable
489  if (phys_energy) then
490  e_c_ = var_set_fluxvar("e_c", "p_c")
491  iw_e = e_c_
492  else
493  e_c_ = -1
494  end if
495 
496  ! ambipolar sts assumes mag and energy charges are continuous
497  allocate(mag(ndir))
498  mag(:) = var_set_bfield(ndir)
499 
500  if (twofl_glm) then
501  psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
502  else
503  psi_ = -1
504  end if
505 
506  ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
507  tweight_c_ = -1
508  if(twofl_trac) then
509  tcoff_c_ = var_set_wextra()
510  iw_tcoff = tcoff_c_
511  if(twofl_trac_type > 2) then
512  tweight_c_ = var_set_wextra()
513  endif
514  else
515  tcoff_c_ = -1
516  end if
517 
518  !now allocate neutrals
519 
520  ! TODO so far number_species is only used to treat them differently
521  ! in the solvers (different cbounds)
522  if (twofl_cbounds_species) then
523  stop_indices(1)=nwflux
524  start_indices(2)=nwflux+1
525  endif
526 
527  ! Determine flux variables
528  rho_n_ = var_set_fluxvar("rho_n", "rho_n")
529  allocate(mom_n(ndir))
530  do idir=1,ndir
531  mom_n(idir) = var_set_fluxvar("m_n","v_n",idir)
532  enddo
533  if (phys_energy) then
534  e_n_ = var_set_fluxvar("e_n", "p_n")
535  else
536  e_n_ = -1
537  end if
538 
539  tweight_n_ = -1
540  if(twofl_trac) then
541  tcoff_n_ = var_set_wextra()
542  if(twofl_trac_type > 2) then
543  tweight_n_ = var_set_wextra()
544  endif
545  else
546  tcoff_n_ = -1
547  end if
548 
549  stop_indices(number_species)=nwflux
550 
551  ! set indices of equi vars and update number_equi_vars
552  number_equi_vars = 0
553  if(has_equi_rho_n0) then
556  endif
557  if(has_equi_pe_n0) then
560  phys_equi_pe=.true.
561  endif
562  if(has_equi_rho_c0) then
565  iw_equi_rho = equi_rho_c0_
566  endif
567  if(has_equi_pe_c0) then
570  iw_equi_p = equi_pe_c0_
571  phys_equi_pe=.true.
572  endif
573 
574  ! set number of variables which need update ghostcells
575  nwgc=nwflux
576 
577  ! determine number of stagger variables
578  nws=ndim
579 
580  ! Check whether custom flux types have been defined
581  if (.not. allocated(flux_type)) then
582  allocate(flux_type(ndir, nw))
583  flux_type = flux_default
584  else if (any(shape(flux_type) /= [ndir, nw])) then
585  call mpistop("phys_check error: flux_type has wrong shape")
586  end if
587 
588  if(ndim>1) then
589  if(twofl_glm) then
590  flux_type(:,psi_)=flux_special
591  do idir=1,ndir
592  flux_type(idir,mag(idir))=flux_special
593  end do
594  else
595  do idir=1,ndir
596  flux_type(idir,mag(idir))=flux_tvdlf
597  end do
598  end if
599  end if
600 
601  phys_get_dt => twofl_get_dt
602  phys_get_cmax => twofl_get_cmax
603  phys_get_a2max => twofl_get_a2max
604  !phys_get_tcutoff => twofl_get_tcutoff_c
605  if(twofl_cbounds_species) then
606  if (mype .eq. 0) print*, "Using different cbounds for each species nspecies = ", number_species
607  phys_get_cbounds => twofl_get_cbounds_species
608  phys_get_h_speed => twofl_get_h_speed_species
609  else
610  if (mype .eq. 0) print*, "Using same cbounds for all species"
611  phys_get_cbounds => twofl_get_cbounds_one
612  phys_get_h_speed => twofl_get_h_speed_one
613  endif
614  phys_get_flux => twofl_get_flux
615  phys_add_source_geom => twofl_add_source_geom
616  phys_add_source => twofl_add_source
617  phys_to_conserved => twofl_to_conserved
618  phys_to_primitive => twofl_to_primitive
619  phys_check_params => twofl_check_params
620  phys_check_w => twofl_check_w
621  phys_write_info => twofl_write_info
622  phys_handle_small_values => twofl_handle_small_values
623  !set equilibrium variables for the new grid
624  if(number_equi_vars>0) then
625  phys_set_equi_vars => set_equi_vars_grid
626  endif
627  ! convert_type is not known here, so associate the corresp. subroutine in check_params
628  if(type_divb==divb_glm) then
629  phys_modify_wlr => twofl_modify_wlr
630  end if
631 
632  ! if using ct stagger grid, boundary divb=0 is not done here
633  if(stagger_grid) then
634  phys_get_ct_velocity => twofl_get_ct_velocity
635  phys_update_faces => twofl_update_faces
636  phys_face_to_center => twofl_face_to_center
637  phys_modify_wlr => twofl_modify_wlr
638  else if(ndim>1) then
639  phys_boundary_adjust => twofl_boundary_adjust
640  end if
641 
642  {^nooned
643  ! clean initial divb
644  if(clean_initial_divb) phys_clean_divb => twofl_clean_divb_multigrid
645  }
646 
647  ! Whether diagonal ghost cells are required for the physics
648  if(type_divb < divb_linde) phys_req_diagonal = .false.
649 
650  ! derive units from basic units
651  call twofl_physical_units()
652 
653  if(.not. phys_energy .and. (twofl_thermal_conduction_c&
654  .or. twofl_thermal_conduction_n)) then
655  call mpistop("thermal conduction needs twofl_energy=T")
656  end if
657 
658  ! initialize thermal conduction module
660  .or. twofl_thermal_conduction_n) then
661  phys_req_diagonal = .true.
662  call sts_init()
664  endif
666  allocate(tc_fl_c)
667  if(has_equi_pe_c0 .and. has_equi_rho_c0) then
668  tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c_with_equi
669  if(phys_internal_e) then
670  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c_with_equi
671  else
672  if(twofl_eq_energy == eq_energy_ki) then
673  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c_with_equi
674  else
675  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c_with_equi
676  endif
677  endif
678  if(twofl_equi_thermal_c) then
679  tc_fl_c%has_equi = .true.
680  tc_fl_c%get_temperature_equi => twofl_get_temperature_c_equi
681  tc_fl_c%get_rho_equi => twofl_get_rho_c_equi
682  else
683  tc_fl_c%has_equi = .false.
684  endif
685  else
686  if(phys_internal_e) then
687  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eint_c
688  else
689  if(twofl_eq_energy == eq_energy_ki) then
690  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_eki_c
691  else
692  tc_fl_c%get_temperature_from_conserved => twofl_get_temperature_from_etot_c
693  endif
694  endif
695  tc_fl_c%get_temperature_from_eint => twofl_get_temperature_from_eint_c
696  endif
697  if(use_twofl_tc_c .eq. mhd_tc) then
700  else if(use_twofl_tc_c .eq. hd_tc) then
703  endif
704  if(.not. phys_internal_e) then
706  endif
708  tc_fl_c%get_rho => get_rhoc_tot
709  tc_fl_c%e_ = e_c_
710  tc_fl_c%Tcoff_ = tcoff_c_
711  end if
713  allocate(tc_fl_n)
715  if(has_equi_pe_n0 .and. has_equi_rho_n0) then
716  tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n_with_equi
717  if(twofl_equi_thermal_n) then
718  tc_fl_n%has_equi = .true.
719  tc_fl_n%get_temperature_equi => twofl_get_temperature_n_equi
720  tc_fl_n%get_rho_equi => twofl_get_rho_n_equi
721  else
722  tc_fl_n%has_equi = .false.
723  endif
724  else
725  tc_fl_n%get_temperature_from_eint => twofl_get_temperature_from_eint_n
726  endif
727  if(phys_internal_e) then
728  if(has_equi_pe_n0 .and. has_equi_rho_n0) then
729  tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n_with_equi
730  else
731  tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_eint_n
732  endif
734  else
735  if(has_equi_pe_n0 .and. has_equi_rho_n0) then
736  tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n_with_equi
737  else
738  tc_fl_n%get_temperature_from_conserved => twofl_get_temperature_from_etot_n
739  endif
742  endif
744  tc_fl_n%get_rho => get_rhon_tot
745  tc_fl_n%e_ = e_n_
746  tc_fl_n%Tcoff_ = tcoff_n_
747  end if
748 
749 
750  if(.not. phys_energy .and. (twofl_radiative_cooling_c&
751  .or. twofl_radiative_cooling_n)) then
752  call mpistop("radiative cooling needs twofl_energy=T")
753  end if
754 
755  if(twofl_equi_thermal .and. (.not. has_equi_pe_c0 .or. .not. has_equi_pe_n0)) then
756  call mpistop("twofl_equi_thermal=T has_equi_pe_n0 and has _equi_pe_c0=T")
757  endif
758 
759  ! initialize thermal conduction module
761  .or. twofl_radiative_cooling_n) then
762  ! Initialize radiative cooling module
763  call radiative_cooling_init_params(twofl_gamma,he_abundance)
765  allocate(rc_fl_c)
767  rc_fl_c%get_rho => get_rhoc_tot
768  rc_fl_c%get_pthermal => twofl_get_pthermal_c
769  rc_fl_c%get_var_Rfactor => rfactor_c
770  rc_fl_c%e_ = e_c_
771  rc_fl_c%Tcoff_ = tcoff_c_
773  rc_fl_c%has_equi = .true.
774  rc_fl_c%get_rho_equi => twofl_get_rho_c_equi
775  rc_fl_c%get_pthermal_equi => twofl_get_pe_c_equi
776  else
777  rc_fl_c%has_equi = .false.
778  end if
779  end if
780  end if
781  allocate(te_fl_c)
782  te_fl_c%get_rho=> get_rhoc_tot
783  te_fl_c%get_pthermal=> twofl_get_pthermal_c
784  te_fl_c%get_var_Rfactor => rfactor_c
785 {^ifthreed
786  phys_te_images => twofl_te_images
787 }
788 
789  ! Initialize viscosity module
790  !!TODO
791  !if (twofl_viscosity) call viscosity_init(phys_wider_stencil,phys_req_diagonal)
792 
793  ! Initialize gravity module
794  if(twofl_gravity) then
795  ! call gravity_init()
797  end if
798 
799  ! Initialize particles module
800  ! For Hall, we need one more reconstructed layer since currents are computed
801  ! in getflux: assuming one additional ghost layer (two for FOURTHORDER) was
802  ! added in nghostcells.
803  if (twofl_hall) then
804  phys_req_diagonal = .true.
805  if (twofl_4th_order) then
806  phys_wider_stencil = 2
807  else
808  phys_wider_stencil = 1
809  end if
810  end if
811 
812  if(twofl_hyperdiffusivity) then
813  allocate(c_shk(1:nwflux))
814  allocate(c_hyp(1:nwflux))
816  end if
817 
818  end subroutine twofl_phys_init
819 
820 {^ifthreed
821  subroutine twofl_te_images
824 
825  select case(convert_type)
826  case('EIvtiCCmpi','EIvtuCCmpi')
828  case('ESvtiCCmpi','ESvtuCCmpi')
830  case('SIvtiCCmpi','SIvtuCCmpi')
832  case('WIvtiCCmpi','WIvtuCCmpi')
834  case default
835  call mpistop("Error in synthesize emission: Unknown convert_type")
836  end select
837  end subroutine twofl_te_images
838 }
839 
840  ! wrappers for STS functions in thermal_conductivity module
841  ! which take as argument the tc_fluid (defined in the physics module)
842  subroutine twofl_sts_set_source_tc_c_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
844  use mod_fix_conserve
846  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
847  double precision, intent(in) :: x(ixI^S,1:ndim)
848  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
849  double precision, intent(in) :: my_dt
850  logical, intent(in) :: fix_conserve_at_step
851  call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_c)
852  end subroutine twofl_sts_set_source_tc_c_mhd
853 
854  subroutine twofl_sts_set_source_tc_c_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
856  use mod_fix_conserve
858  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
859  double precision, intent(in) :: x(ixI^S,1:ndim)
860  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
861  double precision, intent(in) :: my_dt
862  logical, intent(in) :: fix_conserve_at_step
863  call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_c)
864  end subroutine twofl_sts_set_source_tc_c_hd
865 
866  function twofl_get_tc_dt_mhd_c(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
867  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
868  !where tc_k_para_i=tc_k_para*B_i**2/B**2
869  !and T=p/rho
872 
873  integer, intent(in) :: ixi^l, ixo^l
874  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
875  double precision, intent(in) :: w(ixi^s,1:nw)
876  double precision :: dtnew
877 
878  dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl_c)
879  end function twofl_get_tc_dt_mhd_c
880 
881  function twofl_get_tc_dt_hd_c(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
882  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
883  !where tc_k_para_i=tc_k_para*B_i**2/B**2
884  !and T=p/rho
887 
888  integer, intent(in) :: ixi^l, ixo^l
889  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
890  double precision, intent(in) :: w(ixi^s,1:nw)
891  double precision :: dtnew
892 
893  dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl_c)
894  end function twofl_get_tc_dt_hd_c
895 
896  subroutine twofl_tc_handle_small_e_c(w, x, ixI^L, ixO^L, step)
898  use mod_small_values
899 
900  integer, intent(in) :: ixI^L,ixO^L
901  double precision, intent(inout) :: w(ixI^S,1:nw)
902  double precision, intent(in) :: x(ixI^S,1:ndim)
903  integer, intent(in) :: step
904 
905  character(len=140) :: error_msg
906 
907  write(error_msg,"(a,i3)") "Charges thermal conduction step ", step
908  call twofl_handle_small_ei_c(w,x,ixi^l,ixo^l,e_c_,error_msg)
909  end subroutine twofl_tc_handle_small_e_c
910 
911  subroutine twofl_sts_set_source_tc_n_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
913  use mod_fix_conserve
915  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
916  double precision, intent(in) :: x(ixI^S,1:ndim)
917  double precision, intent(inout) :: wres(ixI^S,1:nw), w(ixI^S,1:nw)
918  double precision, intent(in) :: my_dt
919  logical, intent(in) :: fix_conserve_at_step
920  call sts_set_source_tc_hd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl_n)
921  end subroutine twofl_sts_set_source_tc_n_hd
922 
923  subroutine twofl_tc_handle_small_e_n(w, x, ixI^L, ixO^L, step)
925 
926  integer, intent(in) :: ixI^L,ixO^L
927  double precision, intent(inout) :: w(ixI^S,1:nw)
928  double precision, intent(in) :: x(ixI^S,1:ndim)
929  integer, intent(in) :: step
930 
931  character(len=140) :: error_msg
932 
933  write(error_msg,"(a,i3)") "Neutral thermal conduction step ", step
934  call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,error_msg)
935  end subroutine twofl_tc_handle_small_e_n
936 
937  function twofl_get_tc_dt_hd_n(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
938  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
939  !where tc_k_para_i=tc_k_para*B_i**2/B**2
940  !and T=p/rho
943 
944  integer, intent(in) :: ixi^l, ixo^l
945  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
946  double precision, intent(in) :: w(ixi^s,1:nw)
947  double precision :: dtnew
948 
949  dtnew=get_tc_dt_hd(w,ixi^l,ixo^l,dx^d,x,tc_fl_n)
950  end function twofl_get_tc_dt_hd_n
951 
952  subroutine tc_n_params_read_hd(fl)
954  use mod_global_parameters, only: unitpar
955  type(tc_fluid), intent(inout) :: fl
956  integer :: n
957  logical :: tc_saturate=.false.
958  double precision :: tc_k_para=0d0
959 
960  namelist /tc_n_list/ tc_saturate, tc_k_para
961 
962  do n = 1, size(par_files)
963  open(unitpar, file=trim(par_files(n)), status="old")
964  read(unitpar, tc_n_list, end=111)
965 111 close(unitpar)
966  end do
967  fl%tc_saturate = tc_saturate
968  fl%tc_k_para = tc_k_para
969 
970  end subroutine tc_n_params_read_hd
971 
972  subroutine rc_params_read_n(fl)
974  use mod_constants, only: bigdouble
975  type(rc_fluid), intent(inout) :: fl
976  integer :: n
977  ! list parameters
978  integer :: ncool = 4000
979  double precision :: cfrac=0.1d0
980 
981  !> Name of cooling curve
982  character(len=std_len) :: coolcurve='JCorona'
983 
984  !> Name of cooling method
985  character(len=std_len) :: coolmethod='exact'
986 
987  !> Fixed temperature not lower than tlow
988  logical :: Tfix=.false.
989 
990  !> Lower limit of temperature
991  double precision :: tlow=bigdouble
992 
993  !> Add cooling source in a split way (.true.) or un-split way (.false.)
994  logical :: rc_split=.false.
995 
996  namelist /rc_list_n/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
997 
998  do n = 1, size(par_files)
999  open(unitpar, file=trim(par_files(n)), status="old")
1000  read(unitpar, rc_list_n, end=111)
1001 111 close(unitpar)
1002  end do
1003 
1004  fl%ncool=ncool
1005  fl%coolcurve=coolcurve
1006  fl%coolmethod=coolmethod
1007  fl%tlow=tlow
1008  fl%Tfix=tfix
1009  fl%rc_split=rc_split
1010  fl%cfrac=cfrac
1011  end subroutine rc_params_read_n
1012 
1013  !end wrappers
1014 
1015  ! fill in tc_fluid fields from namelist
1016  subroutine tc_c_params_read_mhd(fl)
1018  type(tc_fluid), intent(inout) :: fl
1019 
1020  integer :: n
1021 
1022  ! list parameters
1023  logical :: tc_perpendicular=.false.
1024  logical :: tc_saturate=.false.
1025  double precision :: tc_k_para=0d0
1026  double precision :: tc_k_perp=0d0
1027  character(len=std_len) :: tc_slope_limiter="MC"
1028 
1029  namelist /tc_c_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
1030  do n = 1, size(par_files)
1031  open(unitpar, file=trim(par_files(n)), status="old")
1032  read(unitpar, tc_c_list, end=111)
1033 111 close(unitpar)
1034  end do
1035 
1036  fl%tc_perpendicular = tc_perpendicular
1037  fl%tc_saturate = tc_saturate
1038  fl%tc_k_para = tc_k_para
1039  fl%tc_k_perp = tc_k_perp
1040  select case(tc_slope_limiter)
1041  case ('no','none')
1042  fl%tc_slope_limiter = 0
1043  case ('MC')
1044  ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
1045  fl%tc_slope_limiter = 1
1046  case('minmod')
1047  ! minmod limiter
1048  fl%tc_slope_limiter = 2
1049  case ('superbee')
1050  ! Roes superbee limiter (eq.3.51i)
1051  fl%tc_slope_limiter = 3
1052  case ('koren')
1053  ! Barry Koren Right variant
1054  fl%tc_slope_limiter = 4
1055  case default
1056  call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
1057  end select
1058  end subroutine tc_c_params_read_mhd
1059 
1060  subroutine tc_c_params_read_hd(fl)
1062  use mod_global_parameters, only: unitpar
1063  type(tc_fluid), intent(inout) :: fl
1064  integer :: n
1065  logical :: tc_saturate=.false.
1066  double precision :: tc_k_para=0d0
1067 
1068  namelist /tc_c_list/ tc_saturate, tc_k_para
1069 
1070  do n = 1, size(par_files)
1071  open(unitpar, file=trim(par_files(n)), status="old")
1072  read(unitpar, tc_c_list, end=111)
1073 111 close(unitpar)
1074  end do
1075  fl%tc_saturate = tc_saturate
1076  fl%tc_k_para = tc_k_para
1077 
1078  end subroutine tc_c_params_read_hd
1079 
1080 !! end th cond
1081 
1082 !!rad cool
1083  subroutine rc_params_read_c(fl)
1085  use mod_constants, only: bigdouble
1086  type(rc_fluid), intent(inout) :: fl
1087  integer :: n
1088  ! list parameters
1089  integer :: ncool = 4000
1090  double precision :: cfrac=0.1d0
1091 
1092  !> Name of cooling curve
1093  character(len=std_len) :: coolcurve='JCcorona'
1094 
1095  !> Name of cooling method
1096  character(len=std_len) :: coolmethod='exact'
1097 
1098  !> Fixed temperature not lower than tlow
1099  logical :: Tfix=.false.
1100 
1101  !> Lower limit of temperature
1102  double precision :: tlow=bigdouble
1103 
1104  !> Add cooling source in a split way (.true.) or un-split way (.false.)
1105  logical :: rc_split=.false.
1106 
1107 
1108  namelist /rc_list_c/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split
1109 
1110  do n = 1, size(par_files)
1111  open(unitpar, file=trim(par_files(n)), status="old")
1112  read(unitpar, rc_list_c, end=111)
1113 111 close(unitpar)
1114  end do
1115 
1116  fl%ncool=ncool
1117  fl%coolcurve=coolcurve
1118  fl%coolmethod=coolmethod
1119  fl%tlow=tlow
1120  fl%Tfix=tfix
1121  fl%rc_split=rc_split
1122  fl%cfrac=cfrac
1123  end subroutine rc_params_read_c
1124 
1125 !! end rad cool
1126 
1127  !> sets the equilibrium variables
1128  subroutine set_equi_vars_grid_faces(igrid,x,ixI^L,ixO^L)
1131  use mod_usr_methods
1132  integer, intent(in) :: igrid, ixI^L, ixO^L
1133  double precision, intent(in) :: x(ixI^S,1:ndim)
1134 
1135  double precision :: delx(ixI^S,1:ndim)
1136  double precision :: xC(ixI^S,1:ndim),xshift^D
1137  integer :: idims, ixC^L, hxO^L, ix, idims2
1138 
1139  if(slab_uniform)then
1140  ^d&delx(ixi^s,^d)=rnode(rpdx^d_,igrid)\
1141  else
1142  ! for all non-cartesian and stretched cartesian coordinates
1143  delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
1144  endif
1145 
1146 
1147  do idims=1,ndim
1148  hxo^l=ixo^l-kr(idims,^d);
1149  if(stagger_grid) then
1150  ! ct needs all transverse cells
1151  ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
1152  else
1153  ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
1154  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
1155  end if
1156  ! always xshift=0 or 1/2
1157  xshift^d=half*(one-kr(^d,idims));
1158  do idims2=1,ndim
1159  select case(idims2)
1160  {case(^d)
1161  do ix = ixc^lim^d
1162  ! xshift=half: this is the cell center coordinate
1163  ! xshift=0: this is the cell edge i+1/2 coordinate
1164  xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
1165  end do\}
1166  end select
1167  end do
1168  call usr_set_equi_vars(ixi^l,ixc^l,xc,ps(igrid)%equi_vars(ixi^s,1:number_equi_vars,idims))
1169  end do
1170 
1171  end subroutine set_equi_vars_grid_faces
1172 
1173  !> sets the equilibrium variables
1174  subroutine set_equi_vars_grid(igrid)
1176  use mod_usr_methods
1177 
1178  integer, intent(in) :: igrid
1179 
1180  !values at the center
1181  call usr_set_equi_vars(ixg^ll,ixg^ll,ps(igrid)%x,ps(igrid)%equi_vars(ixg^t,1:number_equi_vars,0))
1182 
1183  !values at the interfaces
1184  call set_equi_vars_grid_faces(igrid,ps(igrid)%x,ixg^ll,ixm^ll)
1185 
1186  end subroutine set_equi_vars_grid
1187 
1188  ! w, wnew conserved
1189  function convert_vars_splitting(ixI^L,ixO^L, w, x, nwc) result(wnew)
1191  integer, intent(in) :: ixi^l,ixo^l, nwc
1192  double precision, intent(in) :: w(ixi^s, 1:nw)
1193  double precision, intent(in) :: x(ixi^s,1:ndim)
1194  double precision :: wnew(ixo^s, 1:nwc)
1195  double precision :: rho(ixi^s)
1196 
1197  call get_rhon_tot(w,x,ixi^l,ixo^l,rho(ixi^s))
1198  wnew(ixo^s,rho_n_) = rho(ixo^s)
1199  wnew(ixo^s,mom_n(:)) = w(ixo^s,mom_n(:))
1200  call get_rhoc_tot(w,x,ixi^l,ixo^l,rho(ixi^s))
1201  wnew(ixo^s,rho_c_) = rho(ixo^s)
1202  wnew(ixo^s,mom_c(:)) = w(ixo^s,mom_c(:))
1203 
1204  if (b0field) then
1205  ! add background magnetic field B0 to B
1206  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))+block%B0(ixo^s,:,0)
1207  else
1208  wnew(ixo^s,mag(:))=w(ixo^s,mag(:))
1209  end if
1210 
1211  if(phys_energy) then
1212  wnew(ixo^s,e_n_) = w(ixo^s,e_n_)
1213  if(has_equi_pe_n0) then
1214  wnew(ixo^s,e_n_) = wnew(ixo^s,e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,0)* inv_gamma_1
1215  endif
1216  wnew(ixo^s,e_c_) = w(ixo^s,e_c_)
1217  if(has_equi_pe_c0) then
1218  wnew(ixo^s,e_c_) = wnew(ixo^s,e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,0)* inv_gamma_1
1219  endif
1220  if(b0field .and. phys_total_energy) then
1221  wnew(ixo^s,e_c_)=wnew(ixo^s,e_c_)+0.5d0*sum(block%B0(ixo^s,:,0)**2,dim=ndim+1) &
1222  + sum(w(ixo^s,mag(:))*block%B0(ixo^s,:,0),dim=ndim+1)
1223  endif
1224  endif
1225 
1226  end function convert_vars_splitting
1227 
1228  !> copied from mod_gravity
1229  subroutine grav_params_read(files)
1230  use mod_global_parameters, only: unitpar
1231  character(len=*), intent(in) :: files(:)
1232  integer :: n
1233 
1234  namelist /grav_list/ grav_split
1235 
1236  do n = 1, size(files)
1237  open(unitpar, file=trim(files(n)), status="old")
1238  read(unitpar, grav_list, end=111)
1239 111 close(unitpar)
1240  end do
1241 
1242  end subroutine grav_params_read
1243 
1246  use mod_convert, only: add_convert_method
1247  integer :: ii
1248  do ii = 1,ndim
1249  if(ii==1) then
1250  call add_convert_method(dump_hyperdiffusivity_coef_x, nw, cons_wnames(1:nw), "hyper_x")
1251  elseif(ii==2) then
1252  call add_convert_method(dump_hyperdiffusivity_coef_y, nw, cons_wnames(1:nw), "hyper_y")
1253  else
1254  call add_convert_method(dump_hyperdiffusivity_coef_z, nw, cons_wnames(1:nw), "hyper_z")
1255  endif
1256  enddo
1257  end subroutine associate_dump_hyper
1258 
1261  use mod_usr_methods
1262  use mod_convert, only: add_convert_method
1263 
1264  ! after user parameter setting
1265  gamma_1=twofl_gamma-1.d0
1266  if (.not. phys_energy) then
1267  if (twofl_gamma <= 0.0d0) call mpistop ("Error: twofl_gamma <= 0")
1268  if (twofl_adiab < 0.0d0) call mpistop ("Error: twofl_adiab < 0")
1270  else
1271  if (twofl_gamma <= 0.0d0 .or. twofl_gamma == 1.0d0) &
1272  call mpistop ("Error: twofl_gamma <= 0 or twofl_gamma == 1")
1273  inv_gamma_1=1.d0/gamma_1
1274  small_e = small_pressure * inv_gamma_1
1275  end if
1276 
1277  ! this has to be done here as use_imex_scheme is not set in init subroutine,
1278  ! but here it is
1279  if(use_imex_scheme) then
1280  if(has_collisions()) then
1281  ! implicit collisional terms update
1282  phys_implicit_update => twofl_implicit_coll_terms_update
1283  phys_evaluate_implicit => twofl_evaluate_implicit
1284  if(mype .eq. 1) then
1285  print*, "IMPLICIT UPDATE with calc_mult_factor", twofl_implicit_calc_mult_method
1286  endif
1287  if(twofl_implicit_calc_mult_method == 1) then
1289  else
1290  calc_mult_factor => calc_mult_factor2
1291  endif
1292  endif
1293  else
1294  ! check dtcoll par for explicit implementation of the coll. terms
1295  if(dtcollpar .le. 0d0 .or. dtcollpar .ge. 1d0) then
1296  if (mype .eq. 0) print*, "Explicit update of coll terms requires 0<dtcollpar<1, dtcollpar set to 0.8."
1297  dtcollpar = 0.8
1298  endif
1299 
1300  endif
1301 ! if(H_ion_fr == 0d0 .and. He_ion_fr == 0d0) then
1302 ! call mpistop("H_ion_fr or He_ion_fr must be > 0 or use hd module")
1303 ! endif
1304 ! if(H_ion_fr == 1d0 .and. He_ion_fr == 1d0) then
1305 ! call mpistop("H_ion_fr or He_ion_fr must be < 1 or use mhd module")
1306 ! endif
1307  if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
1308  call mpistop("usr_set_equi_vars has to be implemented in the user file")
1309  endif
1310  if(convert .or. autoconvert) then
1311  if(convert_type .eq. 'dat_generic_mpi') then
1312  if(twofl_dump_full_vars) then
1313  if(mype .eq. 0) print*, " add conversion method: split -> full "
1314  call add_convert_method(convert_vars_splitting, nw, cons_wnames, "new")
1315  endif
1316  if(twofl_dump_coll_terms) then
1317  if(mype .eq. 0) print*, " add conversion method: dump coll terms "
1318  call add_convert_method(dump_coll_terms, 3, (/"alpha ", "gamma_rec", "gamma_ion"/), "_coll")
1319  endif
1321  if(mype .eq. 0) print*, " add conversion method: dump hyperdiffusivity coeff. "
1322  call associate_dump_hyper()
1323  endif
1324  endif
1325  endif
1326  end subroutine twofl_check_params
1327 
1330  double precision :: mp,kB,miu0,c_lightspeed
1331  !double precision :: a,b,c,d
1332  double precision :: a,b
1333  ! Derive scaling units
1334  if(si_unit) then
1335  mp=mp_si
1336  kb=kb_si
1337  miu0=miu0_si
1338  c_lightspeed=c_si
1339  else
1340  mp=mp_cgs
1341  kb=kb_cgs
1342  miu0=4.d0*dpi
1343  c_lightspeed=const_c
1344  end if
1345 
1346  a=1d0
1347  b=1d0
1348  rc=2d0
1349  rn=1d0
1350 
1351  !now the unit choice:
1352  !unit 1 from number density or density -> mH
1353  !unit 2 from
1354 
1355  if(unit_density/=1.d0) then
1357  else
1358  ! unit of numberdensity is independent by default
1360  end if
1361  if(unit_velocity/=1.d0) then
1365  else if(unit_pressure/=1.d0) then
1369  else if(unit_magneticfield/=1.d0) then
1373  else if(unit_temperature/=1.d0) then
1377  end if
1378  if(unit_time/=1.d0) then
1380  else
1381  ! unit of length is independent by default
1383  end if
1384  ! Additional units needed for the particles
1385  c_norm=c_lightspeed/unit_velocity
1387  if (.not. si_unit) unit_charge = unit_charge*const_c
1389  end subroutine twofl_physical_units
1390 
1391  subroutine twofl_check_w(primitive,ixI^L,ixO^L,w,flag)
1393 
1394  logical, intent(in) :: primitive
1395  integer, intent(in) :: ixI^L, ixO^L
1396  double precision, intent(in) :: w(ixI^S,nw)
1397  double precision :: tmp(ixI^S)
1398  logical, intent(inout) :: flag(ixI^S,1:nw)
1399 
1400  flag=.false.
1401 
1402  if(has_equi_rho_n0) then
1403  tmp(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0)
1404  else
1405  tmp(ixo^s) = w(ixo^s,rho_n_)
1406  endif
1407  where(tmp(ixo^s) < small_density) flag(ixo^s,rho_n_) = .true.
1408  if(has_equi_rho_c0) then
1409  tmp(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0)
1410  else
1411  tmp(ixo^s) = w(ixo^s,rho_c_)
1412  endif
1413  where(tmp(ixo^s) < small_density) flag(ixo^s,rho_c_) = .true.
1414  if(phys_energy) then
1415  if(primitive) then
1416  tmp(ixo^s) = w(ixo^s,e_n_)
1417  if(has_equi_pe_n0) then
1418  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)
1419  endif
1420  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_n_) = .true.
1421  tmp(ixo^s) = w(ixo^s,e_c_)
1422  if(has_equi_pe_c0) then
1423  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)
1424  endif
1425  where(tmp(ixo^s) < small_pressure) flag(ixo^s,e_c_) = .true.
1426  else
1427  if(phys_internal_e) then
1428  tmp(ixo^s)=w(ixo^s,e_n_)
1429  if(has_equi_pe_n0) then
1430  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
1431  endif
1432  where(tmp(ixo^s) < small_e) flag(ixo^s,e_n_) = .true.
1433  tmp(ixo^s)=w(ixo^s,e_c_)
1434  if(has_equi_pe_c0) then
1435  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
1436  endif
1437  where(tmp(ixo^s) < small_e) flag(ixo^s,e_c_) = .true.
1438  else
1439  !neutrals
1440  tmp(ixo^s)=w(ixo^s,e_n_)-&
1441  twofl_kin_en_n(w,ixi^l,ixo^l)
1442  if(has_equi_pe_n0) then
1443  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
1444  endif
1445  where(tmp(ixo^s) < small_e) flag(ixo^s,e_n_) = .true.
1446  if(phys_total_energy) then
1447  tmp(ixo^s)=w(ixo^s,e_c_)-&
1448  twofl_kin_en_c(w,ixi^l,ixo^l)-twofl_mag_en(w,ixi^l,ixo^l)
1449  else
1450  tmp(ixo^s)=w(ixo^s,e_c_)-&
1451  twofl_kin_en_c(w,ixi^l,ixo^l)
1452  end if
1453  if(has_equi_pe_c0) then
1454  tmp(ixo^s) = tmp(ixo^s)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
1455  endif
1456  where(tmp(ixo^s) < small_e) flag(ixo^s,e_c_) = .true.
1457  end if
1458  endif
1459  end if
1460 
1461  end subroutine twofl_check_w
1462 
1463  !> Transform primitive variables into conservative ones
1464  subroutine twofl_to_conserved(ixI^L,ixO^L,w,x)
1466  integer, intent(in) :: ixi^l, ixo^l
1467  double precision, intent(inout) :: w(ixi^s, nw)
1468  double precision, intent(in) :: x(ixi^s, 1:ndim)
1469  integer :: idir
1470  double precision :: rhoc(ixi^s)
1471  double precision :: rhon(ixi^s)
1472 
1473  !if (fix_small_values) then
1474  ! call twofl_handle_small_values(.true., w, x, ixI^L, ixO^L, 'twofl_to_conserved')
1475  !end if
1476 
1477  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
1478  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
1479 
1480  ! Calculate total energy from pressure, kinetic and magnetic energy
1481  if(phys_energy) then
1482  if(phys_internal_e) then
1483  w(ixo^s,e_n_)=w(ixo^s,e_n_)*inv_gamma_1
1484  w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1
1485  else
1486  w(ixo^s,e_n_)=w(ixo^s,e_n_)*inv_gamma_1&
1487  +half*sum(w(ixo^s,mom_n(:))**2,dim=ndim+1)*rhon(ixo^s)
1488  if(phys_total_energy) then
1489  w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1&
1490  +half*sum(w(ixo^s,mom_c(:))**2,dim=ndim+1)*rhoc(ixo^s)&
1491  +twofl_mag_en(w, ixi^l, ixo^l)
1492  else
1493  ! kinetic energy + internal energy is evolved
1494  w(ixo^s,e_c_)=w(ixo^s,e_c_)*inv_gamma_1&
1495  +half*sum(w(ixo^s,mom_c(:))**2,dim=ndim+1)*rhoc(ixo^s)
1496  end if
1497  end if
1498  end if
1499 
1500  ! Convert velocity to momentum
1501  do idir = 1, ndir
1502  w(ixo^s, mom_n(idir)) = rhon(ixo^s) * w(ixo^s, mom_n(idir))
1503  w(ixo^s, mom_c(idir)) = rhoc(ixo^s) * w(ixo^s, mom_c(idir))
1504  end do
1505  end subroutine twofl_to_conserved
1506 
1507  !> Transform conservative variables into primitive ones
1508  subroutine twofl_to_primitive(ixI^L,ixO^L,w,x)
1510  integer, intent(in) :: ixi^l, ixo^l
1511  double precision, intent(inout) :: w(ixi^s, nw)
1512  double precision, intent(in) :: x(ixi^s, 1:ndim)
1513  integer :: idir
1514  double precision :: rhoc(ixi^s)
1515  double precision :: rhon(ixi^s)
1516 
1517  if (fix_small_values) then
1518  call twofl_handle_small_values(.false., w, x, ixi^l, ixo^l, 'twofl_to_primitive')
1519  end if
1520 
1521  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
1522  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
1523 
1524  if(phys_energy) then
1525  if(phys_internal_e) then
1526  w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
1527  w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
1528  else
1529  ! neutrals evolved energy = ke + e_int
1530  w(ixo^s,e_n_)=gamma_1*(w(ixo^s,e_n_)&
1531  -twofl_kin_en_n(w,ixi^l,ixo^l))
1532  ! charges
1533  if(phys_total_energy) then
1534  ! evolved energy = ke + e_int + e_mag
1535  w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
1536  -twofl_kin_en_c(w,ixi^l,ixo^l)&
1537  -twofl_mag_en(w,ixi^l,ixo^l))
1538  else
1539  ! evolved energy = ke + e_int
1540  w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
1541  -twofl_kin_en_c(w,ixi^l,ixo^l))
1542  end if
1543  end if
1544  end if
1545 
1546  ! Convert momentum to velocity
1547  do idir = 1, ndir
1548  w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
1549  w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
1550  end do
1551 
1552  end subroutine twofl_to_primitive
1553 
1554 !!USED IN TC
1555  !> Transform internal energy to total energy
1556  subroutine twofl_ei_to_e_c(ixI^L,ixO^L,w,x)
1558  integer, intent(in) :: ixI^L, ixO^L
1559  double precision, intent(inout) :: w(ixI^S, nw)
1560  double precision, intent(in) :: x(ixI^S, 1:ndim)
1561 
1562  ! Calculate total energy from internal, kinetic and magnetic energy
1563  if(twofl_eq_energy == eq_energy_ki) then
1564  w(ixo^s,e_c_)=w(ixo^s,e_c_)&
1565  +twofl_kin_en_c(w,ixi^l,ixo^l)
1566  else
1567  w(ixo^s,e_c_)=w(ixo^s,e_c_)&
1568  +twofl_kin_en_c(w,ixi^l,ixo^l)&
1569  +twofl_mag_en(w,ixi^l,ixo^l)
1570  endif
1571  end subroutine twofl_ei_to_e_c
1572 
1573  !> Transform total energy to internal energy
1574  subroutine twofl_e_to_ei_c(ixI^L,ixO^L,w,x)
1576  integer, intent(in) :: ixI^L, ixO^L
1577  double precision, intent(inout) :: w(ixI^S, nw)
1578  double precision, intent(in) :: x(ixI^S, 1:ndim)
1579 
1580  if(twofl_eq_energy == eq_energy_ki) then
1581  w(ixo^s,e_c_)=w(ixo^s,e_c_)&
1582  -twofl_kin_en_c(w,ixi^l,ixo^l)
1583  else
1584  ! Calculate ei = e - ek - eb
1585  w(ixo^s,e_c_)=w(ixo^s,e_c_)&
1586  -twofl_kin_en_c(w,ixi^l,ixo^l)&
1587  -twofl_mag_en(w,ixi^l,ixo^l)
1588  endif
1589  end subroutine twofl_e_to_ei_c
1590 
1591  !Neutrals
1592  subroutine twofl_ei_to_e_n(ixI^L,ixO^L,w,x)
1594  integer, intent(in) :: ixI^L, ixO^L
1595  double precision, intent(inout) :: w(ixI^S, nw)
1596  double precision, intent(in) :: x(ixI^S, 1:ndim)
1597 
1598  ! Calculate total energy from internal and kinetic energy
1599 
1600  w(ixo^s,e_n_)=w(ixo^s,e_n_)+twofl_kin_en_n(w,ixi^l,ixo^l)
1601 
1602  end subroutine twofl_ei_to_e_n
1603 
1604  !> Transform total energy to internal energy
1605  subroutine twofl_e_to_ei_n(ixI^L,ixO^L,w,x)
1607  integer, intent(in) :: ixI^L, ixO^L
1608  double precision, intent(inout) :: w(ixI^S, nw)
1609  double precision, intent(in) :: x(ixI^S, 1:ndim)
1610 
1611  ! Calculate ei = e - ek
1612  w(ixo^s,e_n_)=w(ixo^s,e_n_)-twofl_kin_en_n(w,ixi^l,ixo^l)
1613 
1614  call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,"e_to_ei_n")
1615  end subroutine twofl_e_to_ei_n
1616 
1617  subroutine twofl_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
1619  use mod_small_values
1620  logical, intent(in) :: primitive
1621  integer, intent(in) :: ixI^L,ixO^L
1622  double precision, intent(inout) :: w(ixI^S,1:nw)
1623  double precision, intent(in) :: x(ixI^S,1:ndim)
1624  character(len=*), intent(in) :: subname
1625 
1626  integer :: idir
1627  logical :: flag(ixI^S,1:nw)
1628  double precision :: tmp2(ixI^S)
1629  double precision :: tmp1(ixI^S)
1630 
1631  call twofl_check_w(primitive, ixi^l, ixo^l, w, flag)
1632 
1633  if(any(flag)) then
1634  select case (small_values_method)
1635  case ("replace")
1636  if(has_equi_rho_c0) then
1637  where(flag(ixo^s,rho_c_)) w(ixo^s,rho_c_) = &
1638  small_density-block%equi_vars(ixo^s,equi_rho_c0_,0)
1639  else
1640  where(flag(ixo^s,rho_c_)) w(ixo^s,rho_c_) = small_density
1641  end if
1642  if(has_equi_rho_n0) then
1643  where(flag(ixo^s,rho_n_)) w(ixo^s,rho_n_) = &
1644  small_density-block%equi_vars(ixo^s,equi_rho_n0_,0)
1645  else
1646  where(flag(ixo^s,rho_n_)) w(ixo^s,rho_n_) = small_density
1647  end if
1648  do idir = 1, ndir
1649  if(small_values_fix_iw(mom_n(idir))) then
1650  where(flag(ixo^s,rho_n_)) w(ixo^s, mom_n(idir)) = 0.0d0
1651  end if
1652  if(small_values_fix_iw(mom_c(idir))) then
1653  where(flag(ixo^s,rho_c_)) w(ixo^s, mom_c(idir)) = 0.0d0
1654  end if
1655  end do
1656 
1657  if(phys_energy) then
1658  if(primitive) then
1659  if(has_equi_pe_n0) then
1660  tmp1(ixo^s) = small_pressure - &
1661  block%equi_vars(ixo^s,equi_pe_n0_,0)
1662  else
1663  tmp1(ixo^s) = small_pressure
1664  end if
1665  if(has_equi_pe_c0) then
1666  tmp2(ixo^s) = small_e - &
1667  block%equi_vars(ixo^s,equi_pe_c0_,0)
1668  else
1669  tmp2(ixo^s) = small_pressure
1670  end if
1671  else
1672  ! conserved
1673  if(has_equi_pe_n0) then
1674  tmp1(ixo^s) = small_e - &
1675  block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
1676  else
1677  tmp1(ixo^s) = small_e
1678  end if
1679  if(has_equi_pe_c0) then
1680  tmp2(ixo^s) = small_e - &
1681  block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
1682  else
1683  tmp2(ixo^s) = small_e
1684  end if
1685  if(phys_internal_e) then
1686  where(flag(ixo^s,e_n_))
1687  w(ixo^s,e_n_)=tmp1(ixo^s)
1688  end where
1689  where(flag(ixo^s,e_c_))
1690  w(ixo^s,e_c_)=tmp2(ixo^s)
1691  end where
1692  else
1693  where(flag(ixo^s,e_n_))
1694  w(ixo^s,e_n_) = tmp1(ixo^s)+&
1695  twofl_kin_en_n(w,ixi^l,ixo^l)
1696  end where
1697  if(phys_total_energy) then
1698  where(flag(ixo^s,e_c_))
1699  w(ixo^s,e_c_) = tmp2(ixo^s)+&
1700  twofl_kin_en_c(w,ixi^l,ixo^l)+&
1701  twofl_mag_en(w,ixi^l,ixo^l)
1702  end where
1703  else
1704  where(flag(ixo^s,e_c_))
1705  w(ixo^s,e_c_) = tmp2(ixo^s)+&
1706  twofl_kin_en_c(w,ixi^l,ixo^l)
1707  end where
1708  end if
1709  end if
1710  end if
1711  end if
1712  case ("average")
1713  call small_values_average(ixi^l, ixo^l, w, x, flag)
1714  case default
1715  if(.not.primitive) then
1716  !convert w to primitive
1717  ! Calculate pressure = (gamma-1) * (e-ek-eb)
1718  if(phys_energy) then
1719  if(phys_internal_e) then
1720  w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
1721  w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
1722  else
1723  w(ixo^s,e_n_)=gamma_1*(w(ixo^s,e_n_)&
1724  -twofl_kin_en_n(w,ixi^l,ixo^l))
1725  if(phys_total_energy) then
1726  w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
1727  -twofl_kin_en_c(w,ixi^l,ixo^l)&
1728  -twofl_mag_en(w,ixi^l,ixo^l))
1729  else
1730  w(ixo^s,e_c_)=gamma_1*(w(ixo^s,e_c_)&
1731  -twofl_kin_en_c(w,ixi^l,ixo^l))
1732 
1733  end if
1734  end if
1735  end if
1736  ! Convert momentum to velocity
1737  if(has_equi_rho_n0) then
1738  tmp1(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0)
1739  else
1740  tmp1(ixo^s) = w(ixo^s,rho_n_)
1741  end if
1742 
1743  if(has_equi_rho_c0) then
1744  tmp2(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0)
1745  else
1746  tmp2(ixo^s) = w(ixo^s,rho_c_)
1747  end if
1748  do idir = 1, ndir
1749  w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/tmp1(ixo^s)
1750  w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/tmp2(ixo^s)
1751  end do
1752  end if
1753  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1754  end select
1755  end if
1756  end subroutine twofl_handle_small_values
1757 
1758  !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
1759  subroutine twofl_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
1761 
1762  integer, intent(in) :: ixI^L, ixO^L, idim
1763  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1764  double precision, intent(inout) :: cmax(ixI^S)
1765  double precision :: vc(ixI^S)
1766  double precision :: cmax2(ixI^S)
1767  double precision :: vn(ixI^S)
1768 
1769  call twofl_get_csound_c_idim(w,x,ixi^l,ixo^l,idim,cmax)
1770  call twofl_get_v_c_idim(w,x,ixi^l,ixo^l,idim,vc)
1771  call twofl_get_v_n_idim(w,x,ixi^l,ixo^l,idim,vn)
1772  call twofl_get_csound_n(w,x,ixi^l,ixo^l,cmax2)
1773  cmax(ixo^s)=max(abs(vn(ixo^s))+cmax2(ixo^s),&
1774  abs(vc(ixo^s))+cmax(ixo^s))
1775 
1776  end subroutine twofl_get_cmax
1777 
1778  subroutine twofl_get_a2max(w,x,ixI^L,ixO^L,a2max)
1780 
1781  integer, intent(in) :: ixI^L, ixO^L
1782  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1783  double precision, intent(inout) :: a2max(ndim)
1784  double precision :: a2(ixI^S,ndim,nw)
1785  integer :: gxO^L,hxO^L,jxO^L,kxO^L,i,j
1786 
1787  a2=zero
1788  do i = 1,ndim
1789  !> 4th order
1790  hxo^l=ixo^l-kr(i,^d);
1791  gxo^l=hxo^l-kr(i,^d);
1792  jxo^l=ixo^l+kr(i,^d);
1793  kxo^l=jxo^l+kr(i,^d);
1794  a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
1795  -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
1796  a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
1797  end do
1798  end subroutine twofl_get_a2max
1799 
1800  ! COPIED from hd/moh_hd_phys
1801  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
1802  subroutine twofl_get_tcutoff_n(ixI^L,ixO^L,w,x,tco_local,Tmax_local)
1804  integer, intent(in) :: ixI^L,ixO^L
1805  double precision, intent(in) :: x(ixI^S,1:ndim),w(ixI^S,1:nw)
1806  double precision, intent(out) :: tco_local, Tmax_local
1807 
1808  double precision, parameter :: delta=0.25d0
1809  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1810  integer :: jxO^L,hxO^L
1811  logical :: lrlt(ixI^S)
1812 
1813  {^ifoned
1814  ! reuse lts as rhon
1815  call get_rhon_tot(w,x,ixi^l,ixi^l,lts)
1816  tmp1(ixi^s)=w(ixi^s,e_n_)-0.5d0*sum(w(ixi^s,mom_n(:))**2,dim=ndim+1)/lts(ixi^s)
1817  te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(twofl_gamma-1.d0)
1818 
1819  tmax_local=maxval(te(ixo^s))
1820 
1821  hxo^l=ixo^l-1;
1822  jxo^l=ixo^l+1;
1823  lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1824  lrlt=.false.
1825  where(lts(ixo^s) > delta)
1826  lrlt(ixo^s)=.true.
1827  end where
1828  tco_local=zero
1829  if(any(lrlt(ixo^s))) then
1830  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1831  end if
1832  }
1833  end subroutine twofl_get_tcutoff_n
1834 
1835  !> get adaptive cutoff temperature for TRAC (Johnston 2019 ApJL, 873, L22)
1836  subroutine twofl_get_tcutoff_c(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
1838  use mod_geometry
1839  integer, intent(in) :: ixI^L,ixO^L
1840  double precision, intent(in) :: x(ixI^S,1:ndim)
1841  double precision, intent(inout) :: w(ixI^S,1:nw)
1842  double precision, intent(out) :: Tco_local,Tmax_local
1843 
1844  double precision, parameter :: trac_delta=0.25d0
1845  double precision :: tmp1(ixI^S),Te(ixI^S),lts(ixI^S)
1846  double precision, dimension(ixI^S,1:ndir) :: bunitvec
1847  double precision, dimension(ixI^S,1:ndim) :: gradT
1848  double precision :: Bdir(ndim)
1849  double precision :: ltr(ixI^S),ltrc,ltrp,altr(ixI^S)
1850  integer :: idims,jxO^L,hxO^L,ixA^D,ixB^D
1851  integer :: jxP^L,hxP^L,ixP^L
1852  logical :: lrlt(ixI^S)
1853 
1854  ! reuse lts as rhoc
1855  call get_rhoc_tot(w,x,ixi^l,ixi^l,lts)
1856  if(phys_internal_e) then
1857  tmp1(ixi^s)=w(ixi^s,e_c_)
1858  else
1859  tmp1(ixi^s)=w(ixi^s,e_c_)-0.5d0*(sum(w(ixi^s,mom_c(:))**2,dim=ndim+1)/&
1860  lts(ixi^s)+sum(w(ixi^s,mag(:))**2,dim=ndim+1))
1861  end if
1862  te(ixi^s)=tmp1(ixi^s)/lts(ixi^s)*(twofl_gamma-1.d0)
1863  tmax_local=maxval(te(ixo^s))
1864 
1865  {^ifoned
1866  select case(twofl_trac_type)
1867  case(0)
1868  !> test case, fixed cutoff temperature
1869  w(ixi^s,tcoff_c_)=2.5d5/unit_temperature
1870  case(1)
1871  hxo^l=ixo^l-1;
1872  jxo^l=ixo^l+1;
1873  lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
1874  lrlt=.false.
1875  where(lts(ixo^s) > trac_delta)
1876  lrlt(ixo^s)=.true.
1877  end where
1878  if(any(lrlt(ixo^s))) then
1879  tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
1880  end if
1881  case(2)
1882  !> iijima et al. 2021, LTRAC method
1883  ltrc=1.5d0
1884  ltrp=2.5d0
1885  ixp^l=ixo^l^ladd1;
1886  hxo^l=ixo^l-1;
1887  jxo^l=ixo^l+1;
1888  hxp^l=ixp^l-1;
1889  jxp^l=ixp^l+1;
1890  lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
1891  ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1892  w(ixo^s,tcoff_c_)=te(ixo^s)*&
1893  (0.25*(ltr(jxo^s)+two*ltr(ixo^s)+ltr(hxo^s)))**0.4d0
1894  case default
1895  call mpistop("twofl_trac_type not allowed for 1D simulation")
1896  end select
1897  }
1898  {^nooned
1899  select case(twofl_trac_type)
1900  case(0)
1901  !> test case, fixed cutoff temperature
1902  w(ixi^s,tcoff_c_)=2.5d5/unit_temperature
1903  case(1,4,6)
1904  ! temperature gradient at cell centers
1905  do idims=1,ndim
1906  call gradient(te,ixi^l,ixo^l,idims,tmp1)
1907  gradt(ixo^s,idims)=tmp1(ixo^s)
1908  end do
1909  ! B vector
1910  if(b0field) then
1911  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
1912  else
1913  bunitvec(ixo^s,:)=w(ixo^s,iw_mag(:))
1914  end if
1915  if(twofl_trac_type .gt. 1) then
1916  ! B direction at cell center
1917  bdir=zero
1918  {do ixa^d=0,1\}
1919  ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
1920  bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
1921  {end do\}
1922  if(sum(bdir(:)**2) .gt. zero) then
1923  bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
1924  end if
1925  block%special_values(3:ndim+2)=bdir(1:ndim)
1926  end if
1927  tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
1928  where(tmp1(ixo^s)/=0.d0)
1929  tmp1(ixo^s)=1.d0/tmp1(ixo^s)
1930  elsewhere
1931  tmp1(ixo^s)=bigdouble
1932  end where
1933  ! b unit vector: magnetic field direction vector
1934  do idims=1,ndim
1935  bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
1936  end do
1937  ! temperature length scale inversed
1938  lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
1939  ! fraction of cells size to temperature length scale
1940  if(slab_uniform) then
1941  lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
1942  else
1943  lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
1944  end if
1945  lrlt=.false.
1946  where(lts(ixo^s) > trac_delta)
1947  lrlt(ixo^s)=.true.
1948  end where
1949  if(any(lrlt(ixo^s))) then
1950  block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
1951  else
1952  block%special_values(1)=zero
1953  end if
1954  block%special_values(2)=tmax_local
1955  case(2)
1956  !> iijima et al. 2021, LTRAC method
1957  ltrc=1.5d0
1958  ltrp=4.d0
1959  ixp^l=ixo^l^ladd1;
1960  ! temperature gradient at cell centers
1961  do idims=1,ndim
1962  call gradient(te,ixi^l,ixp^l,idims,tmp1)
1963  gradt(ixp^s,idims)=tmp1(ixp^s)
1964  end do
1965  ! B vector
1966  if(b0field) then
1967  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))+block%B0(ixp^s,:,0)
1968  else
1969  bunitvec(ixp^s,:)=w(ixp^s,iw_mag(:))
1970  end if
1971  tmp1(ixp^s)=dsqrt(sum(bunitvec(ixp^s,:)**2,dim=ndim+1))
1972  where(tmp1(ixp^s)/=0.d0)
1973  tmp1(ixp^s)=1.d0/tmp1(ixp^s)
1974  elsewhere
1975  tmp1(ixp^s)=bigdouble
1976  end where
1977  ! b unit vector: magnetic field direction vector
1978  do idims=1,ndim
1979  bunitvec(ixp^s,idims)=bunitvec(ixp^s,idims)*tmp1(ixp^s)
1980  end do
1981  ! temperature length scale inversed
1982  lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1983  ! fraction of cells size to temperature length scale
1984  if(slab_uniform) then
1985  lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1986  else
1987  lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1988  end if
1989  ltr(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1990 
1991  altr(ixi^s)=zero
1992  do idims=1,ndim
1993  hxo^l=ixo^l-kr(idims,^d);
1994  jxo^l=ixo^l+kr(idims,^d);
1995  altr(ixo^s)=altr(ixo^s) &
1996  +0.25*(ltr(hxo^s)+two*ltr(ixo^s)+ltr(jxo^s))*bunitvec(ixo^s,idims)**2
1997  w(ixo^s,tcoff_c_)=te(ixo^s)*altr(ixo^s)**(0.4*ltrp)
1998  end do
1999  case(3,5)
2000  !> do nothing here
2001  case default
2002  call mpistop("unknown twofl_trac_type")
2003  end select
2004  }
2005  end subroutine twofl_get_tcutoff_c
2006 
2007  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
2008  subroutine twofl_get_h_speed_one(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2010 
2011  integer, intent(in) :: ixI^L, ixO^L, idim
2012  double precision, intent(in) :: wprim(ixI^S, nw)
2013  double precision, intent(in) :: x(ixI^S,1:ndim)
2014  double precision, intent(out) :: Hspeed(ixI^S,1:number_species)
2015 
2016  double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2017  integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2018 
2019  hspeed=0.d0
2020  ixa^l=ixo^l^ladd1;
2021  do id=1,ndim
2022  call twofl_get_csound_prim(wprim,x,ixi^l,ixa^l,id,tmp)
2023  csound(ixa^s,id)=tmp(ixa^s)
2024  end do
2025  ixcmax^d=ixomax^d;
2026  ixcmin^d=ixomin^d+kr(idim,^d)-1;
2027  jxcmax^d=ixcmax^d+kr(idim,^d);
2028  jxcmin^d=ixcmin^d+kr(idim,^d);
2029  hspeed(ixc^s,1)=0.5d0*abs(&
2030  0.5d0 * (wprim(jxc^s,mom_c(idim))+ wprim(jxc^s,mom_n(idim))) &
2031  +csound(jxc^s,idim)- &
2032  0.5d0 * (wprim(ixc^s,mom_c(idim)) + wprim(ixc^s,mom_n(idim)))&
2033  +csound(ixc^s,idim))
2034 
2035  do id=1,ndim
2036  if(id==idim) cycle
2037  ixamax^d=ixcmax^d+kr(id,^d);
2038  ixamin^d=ixcmin^d+kr(id,^d);
2039  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2040  0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
2041  +csound(ixa^s,id)-&
2042  0.5d0 * (wprim(ixc^s,mom_c(id)) + wprim(ixc^s,mom_n(id)))&
2043  +csound(ixc^s,id)))
2044 
2045 
2046  ixamax^d=ixcmax^d-kr(id,^d);
2047  ixamin^d=ixcmin^d-kr(id,^d);
2048  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2049  0.5d0 * (wprim(ixc^s,mom_c(id)) + wprim(ixc^s,mom_n(id)))&
2050  +csound(ixc^s,id)-&
2051  0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
2052  +csound(ixa^s,id)))
2053 
2054  end do
2055 
2056  do id=1,ndim
2057  if(id==idim) cycle
2058  ixamax^d=jxcmax^d+kr(id,^d);
2059  ixamin^d=jxcmin^d+kr(id,^d);
2060  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2061  0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
2062  +csound(ixa^s,id)-&
2063  0.5d0 * (wprim(jxc^s,mom_c(id)) + wprim(jxc^s,mom_n(id)))&
2064  +csound(jxc^s,id)))
2065  ixamax^d=jxcmax^d-kr(id,^d);
2066  ixamin^d=jxcmin^d-kr(id,^d);
2067  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(&
2068  0.5d0 * (wprim(jxc^s,mom_c(id)) + wprim(jxc^s,mom_n(id)))&
2069  +csound(jxc^s,id)-&
2070  0.5d0 * (wprim(ixa^s,mom_c(id)) + wprim(ixa^s,mom_n(id)))&
2071  +csound(ixa^s,id)))
2072  end do
2073 
2074  end subroutine twofl_get_h_speed_one
2075 
2076  !> get H speed for H-correction to fix the carbuncle problem at grid-aligned shock front
2077  subroutine twofl_get_h_speed_species(wprim,x,ixI^L,ixO^L,idim,Hspeed)
2079 
2080  integer, intent(in) :: ixI^L, ixO^L, idim
2081  double precision, intent(in) :: wprim(ixI^S, nw)
2082  double precision, intent(in) :: x(ixI^S,1:ndim)
2083  double precision, intent(out) :: Hspeed(ixI^S,1:number_species)
2084 
2085  double precision :: csound(ixI^S,ndim),tmp(ixI^S)
2086  integer :: jxC^L, ixC^L, ixA^L, id, ix^D
2087 
2088  hspeed=0.d0
2089  ! charges
2090  ixa^l=ixo^l^ladd1;
2091  do id=1,ndim
2092  call twofl_get_csound_prim_c(wprim,x,ixi^l,ixa^l,id,tmp)
2093  csound(ixa^s,id)=tmp(ixa^s)
2094  end do
2095  ixcmax^d=ixomax^d;
2096  ixcmin^d=ixomin^d+kr(idim,^d)-1;
2097  jxcmax^d=ixcmax^d+kr(idim,^d);
2098  jxcmin^d=ixcmin^d+kr(idim,^d);
2099  hspeed(ixc^s,1)=0.5d0*abs(wprim(jxc^s,mom_c(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom_c(idim))+csound(ixc^s,idim))
2100 
2101  do id=1,ndim
2102  if(id==idim) cycle
2103  ixamax^d=ixcmax^d+kr(id,^d);
2104  ixamin^d=ixcmin^d+kr(id,^d);
2105  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom_c(id))+csound(ixa^s,id)-wprim(ixc^s,mom_c(id))+csound(ixc^s,id)))
2106  ixamax^d=ixcmax^d-kr(id,^d);
2107  ixamin^d=ixcmin^d-kr(id,^d);
2108  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixc^s,mom_c(id))+csound(ixc^s,id)-wprim(ixa^s,mom_c(id))+csound(ixa^s,id)))
2109  end do
2110 
2111  do id=1,ndim
2112  if(id==idim) cycle
2113  ixamax^d=jxcmax^d+kr(id,^d);
2114  ixamin^d=jxcmin^d+kr(id,^d);
2115  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(ixa^s,mom_c(id))+csound(ixa^s,id)-wprim(jxc^s,mom_c(id))+csound(jxc^s,id)))
2116  ixamax^d=jxcmax^d-kr(id,^d);
2117  ixamin^d=jxcmin^d-kr(id,^d);
2118  hspeed(ixc^s,1)=max(hspeed(ixc^s,1),0.5d0*abs(wprim(jxc^s,mom_c(id))+csound(jxc^s,id)-wprim(ixa^s,mom_c(id))+csound(ixa^s,id)))
2119  end do
2120 
2121  ! neutrals
2122  ixa^l=ixo^l^ladd1;
2123  do id=1,ndim
2124  call twofl_get_csound_prim_n(wprim,x,ixi^l,ixa^l,id,tmp)
2125  csound(ixa^s,id)=tmp(ixa^s)
2126  end do
2127  ixcmax^d=ixomax^d;
2128  ixcmin^d=ixomin^d+kr(idim,^d)-1;
2129  jxcmax^d=ixcmax^d+kr(idim,^d);
2130  jxcmin^d=ixcmin^d+kr(idim,^d);
2131  hspeed(ixc^s,2)=0.5d0*abs(wprim(jxc^s,mom_n(idim))+csound(jxc^s,idim)-wprim(ixc^s,mom_n(idim))+csound(ixc^s,idim))
2132 
2133  do id=1,ndim
2134  if(id==idim) cycle
2135  ixamax^d=ixcmax^d+kr(id,^d);
2136  ixamin^d=ixcmin^d+kr(id,^d);
2137  hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,mom_n(id))+csound(ixa^s,id)-wprim(ixc^s,mom_n(id))+csound(ixc^s,id)))
2138  ixamax^d=ixcmax^d-kr(id,^d);
2139  ixamin^d=ixcmin^d-kr(id,^d);
2140  hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixc^s,mom_n(id))+csound(ixc^s,id)-wprim(ixa^s,mom_n(id))+csound(ixa^s,id)))
2141  end do
2142 
2143  do id=1,ndim
2144  if(id==idim) cycle
2145  ixamax^d=jxcmax^d+kr(id,^d);
2146  ixamin^d=jxcmin^d+kr(id,^d);
2147  hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(ixa^s,mom_n(id))+csound(ixa^s,id)-wprim(jxc^s,mom_n(id))+csound(jxc^s,id)))
2148  ixamax^d=jxcmax^d-kr(id,^d);
2149  ixamin^d=jxcmin^d-kr(id,^d);
2150  hspeed(ixc^s,2)=max(hspeed(ixc^s,2),0.5d0*abs(wprim(jxc^s,mom_n(id))+csound(jxc^s,id)-wprim(ixa^s,mom_n(id))+csound(ixa^s,id)))
2151  end do
2152 
2153  end subroutine twofl_get_h_speed_species
2154 
2155  !> Estimating bounds for the minimum and maximum signal velocities
2156  subroutine twofl_get_cbounds_one(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2159 
2160  integer, intent(in) :: ixI^L, ixO^L, idim
2161  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2162  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2163  double precision, intent(in) :: x(ixI^S,1:ndim)
2164  double precision, intent(inout) :: cmax(ixI^S,number_species)
2165  double precision, intent(inout), optional :: cmin(ixI^S,number_species)
2166  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2167 
2168  double precision :: wmean(ixI^S,nw)
2169  double precision :: rhon(ixI^S)
2170  double precision :: rhoc(ixI^S)
2171  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2172  integer :: ix^D
2173 
2174  select case (boundspeed)
2175  case (1)
2176  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2177  ! Methods for Fluid Dynamics" by Toro.
2178  call get_rhoc_tot(wlp,x,ixi^l,ixo^l,rhoc)
2179  call get_rhon_tot(wlp,x,ixi^l,ixo^l,rhon)
2180  tmp1(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2181 
2182  call get_rhoc_tot(wrp,x,ixi^l,ixo^l,rhoc)
2183  call get_rhon_tot(wrp,x,ixi^l,ixo^l,rhon)
2184  tmp2(ixo^s)=sqrt(abs(rhoc(ixo^s) +rhon(ixo^s)))
2185 
2186  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2187  umean(ixo^s)=(0.5*(wlp(ixo^s,mom_n(idim))+wlp(ixo^s,mom_c(idim)))*tmp1(ixo^s) + &
2188  0.5*(wrp(ixo^s,mom_n(idim))+wrp(ixo^s,mom_c(idim)))*tmp2(ixo^s))*tmp3(ixo^s)
2189  call twofl_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2190  call twofl_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2191 
2192  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2193  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*(&
2194  0.5*(wrp(ixo^s,mom_n(idim))+wrp(ixo^s,mom_c(idim)))- &
2195  0.5*(wlp(ixo^s,mom_n(idim))+wlp(ixo^s,mom_c(idim))))**2
2196  dmean(ixo^s)=sqrt(dmean(ixo^s))
2197  if(present(cmin)) then
2198  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2199  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2200  if(h_correction) then
2201  {do ix^db=ixomin^db,ixomax^db\}
2202  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2203  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2204  {end do\}
2205  end if
2206  else
2207  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2208  end if
2209  case (2)
2210  ! typeboundspeed=='cmaxmean'
2211  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2212  call get_rhon_tot(wmean,x,ixi^l,ixo^l,rhon)
2213  tmp2(ixo^s)=wmean(ixo^s,mom_n(idim))/rhon(ixo^s)
2214  call get_rhoc_tot(wmean,x,ixi^l,ixo^l,rhoc)
2215  tmp1(ixo^s)=wmean(ixo^s,mom_c(idim))/rhoc(ixo^s)
2216  call twofl_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
2217  if(present(cmin)) then
2218  cmax(ixo^s,1)=max(max(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) +csoundr(ixo^s),zero)
2219  cmin(ixo^s,1)=min(min(abs(tmp2(ixo^s)), abs(tmp1(ixo^s)) ) -csoundr(ixo^s),zero)
2220  if(h_correction) then
2221  {do ix^db=ixomin^db,ixomax^db\}
2222  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2223  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2224  {end do\}
2225  end if
2226  else
2227  cmax(ixo^s,1)= max(abs(tmp2(ixo^s)),abs(tmp1(ixo^s)))+csoundr(ixo^s)
2228  end if
2229  case (3)
2230  ! Miyoshi 2005 JCP 208, 315 equation (67)
2231  call twofl_get_csound(wlp,x,ixi^l,ixo^l,idim,csoundl)
2232  call twofl_get_csound(wrp,x,ixi^l,ixo^l,idim,csoundr)
2233  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2234  if(present(cmin)) then
2235  cmin(ixo^s,1)=min(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
2236  0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))-csoundl(ixo^s)
2237  cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
2238  0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))+csoundl(ixo^s)
2239  if(h_correction) then
2240  {do ix^db=ixomin^db,ixomax^db\}
2241  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2242  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2243  {end do\}
2244  end if
2245  else
2246  cmax(ixo^s,1)=max(0.5*(wlp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))),&
2247  0.5*(wrp(ixo^s,mom_c(idim))+ wrp(ixo^s,mom_n(idim))))+csoundl(ixo^s)
2248  end if
2249  end select
2250 
2251  end subroutine twofl_get_cbounds_one
2252 
2253  !> Calculate fast magnetosonic wave speed
2254  subroutine twofl_get_csound_prim_c(w,x,ixI^L,ixO^L,idim,csound)
2256 
2257  integer, intent(in) :: ixI^L, ixO^L, idim
2258  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2259  double precision, intent(out):: csound(ixI^S)
2260  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2261  double precision :: inv_rho(ixO^S)
2262  double precision :: rhoc(ixI^S)
2263 
2264  integer :: ix1,ix2
2265 
2266 
2267  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
2268  inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2269 
2270  if(phys_energy) then
2271  call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,csound)
2272  csound(ixo^s)=twofl_gamma*csound(ixo^s)/rhoc(ixo^s)
2273  else
2274  call twofl_get_csound2_adiab_c(w,x,ixi^l,ixo^l,csound)
2275  endif
2276 
2277  ! store |B|^2 in v
2278  b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2279  cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2280  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2281  * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2282  * inv_rho(ixo^s)
2283 
2284  where(avmincs2(ixo^s)<zero)
2285  avmincs2(ixo^s)=zero
2286  end where
2287 
2288  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2289 
2290  if (.not. twofl_hall) then
2291  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2292  else
2293  ! take the Hall velocity into account:
2294  ! most simple estimate, high k limit:
2295  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2296  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2297  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2298  twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2299  end if
2300 
2301  end subroutine twofl_get_csound_prim_c
2302 
2303  !> Calculate fast magnetosonic wave speed
2304  subroutine twofl_get_csound_prim_n(w,x,ixI^L,ixO^L,idim,csound)
2306 
2307  integer, intent(in) :: ixI^L, ixO^L, idim
2308  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2309  double precision, intent(out):: csound(ixI^S)
2310  double precision :: rhon(ixI^S)
2311 
2312  if(phys_energy) then
2313  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2314  call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound)
2315  csound(ixo^s)=twofl_gamma*csound(ixo^s)/rhon(ixo^s)
2316  else
2317  call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound)
2318  endif
2319  csound(ixo^s) = sqrt(csound(ixo^s))
2320 
2321  end subroutine twofl_get_csound_prim_n
2322 
2323  !> Estimating bounds for the minimum and maximum signal velocities
2324  subroutine twofl_get_cbounds_species(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
2327  use mod_variables
2328 
2329  integer, intent(in) :: ixI^L, ixO^L, idim
2330  double precision, intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
2331  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2332  double precision, intent(in) :: x(ixI^S,1:ndim)
2333  double precision, intent(inout) :: cmax(ixI^S,1:number_species)
2334  double precision, intent(inout), optional :: cmin(ixI^S,1:number_species)
2335  double precision, intent(in) :: Hspeed(ixI^S,1:number_species)
2336 
2337  double precision :: wmean(ixI^S,nw)
2338  double precision :: rho(ixI^S)
2339  double precision, dimension(ixI^S) :: umean, dmean, csoundL, csoundR, tmp1,tmp2,tmp3
2340  integer :: ix^D
2341 
2342  select case (boundspeed)
2343  case (1)
2344  ! This implements formula (10.52) from "Riemann Solvers and Numerical
2345  ! Methods for Fluid Dynamics" by Toro.
2346  ! charges
2347  call get_rhoc_tot(wlp,x,ixi^l,ixo^l,rho)
2348  tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2349 
2350  call get_rhoc_tot(wrp,x,ixi^l,ixo^l,rho)
2351  tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2352 
2353  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2354  umean(ixo^s)=(wlp(ixo^s,mom_c(idim))*tmp1(ixo^s)+wrp(ixo^s,mom_c(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2355  call twofl_get_csound_prim_c(wlp,x,ixi^l,ixo^l,idim,csoundl)
2356  call twofl_get_csound_prim_c(wrp,x,ixi^l,ixo^l,idim,csoundr)
2357 
2358 
2359  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2360  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2361  (wrp(ixo^s,mom_c(idim)) - wlp(ixo^s,mom_c(idim)))**2
2362  dmean(ixo^s)=sqrt(dmean(ixo^s))
2363  if(present(cmin)) then
2364  cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
2365  cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
2366  if(h_correction) then
2367  {do ix^db=ixomin^db,ixomax^db\}
2368  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2369  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2370  {end do\}
2371  end if
2372  else
2373  cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
2374  end if
2375 
2376  ! neutrals
2377 
2378  call get_rhon_tot(wlp,x,ixi^l,ixo^l,rho)
2379  tmp1(ixo^s)=sqrt(abs(rho(ixo^s)))
2380 
2381  call get_rhon_tot(wrp,x,ixi^l,ixo^l,rho)
2382  tmp2(ixo^s)=sqrt(abs(rho(ixo^s)))
2383 
2384  tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
2385  umean(ixo^s)=(wlp(ixo^s,mom_n(idim))*tmp1(ixo^s)+wrp(ixo^s,mom_n(idim))*tmp2(ixo^s))*tmp3(ixo^s)
2386  call twofl_get_csound_prim_n(wlp,x,ixi^l,ixo^l,idim,csoundl)
2387  call twofl_get_csound_prim_n(wrp,x,ixi^l,ixo^l,idim,csoundr)
2388 
2389 
2390  dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)**2+tmp2(ixo^s)*csoundr(ixo^s)**2)*tmp3(ixo^s)+&
2391  0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2*&
2392  (wrp(ixo^s,mom_n(idim)) - wlp(ixo^s,mom_n(idim)))**2
2393  dmean(ixo^s)=sqrt(dmean(ixo^s))
2394  if(present(cmin)) then
2395  cmin(ixo^s,2)=umean(ixo^s)-dmean(ixo^s)
2396  cmax(ixo^s,2)=umean(ixo^s)+dmean(ixo^s)
2397  if(h_correction) then
2398  {do ix^db=ixomin^db,ixomax^db\}
2399  cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2400  cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2401  {end do\}
2402  end if
2403  else
2404  cmax(ixo^s,2)=abs(umean(ixo^s))+dmean(ixo^s)
2405  end if
2406 
2407  case (2)
2408  ! typeboundspeed=='cmaxmean'
2409  wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
2410  ! charges
2411 
2412  call get_rhoc_tot(wmean,x,ixi^l,ixo^l,rho)
2413  tmp1(ixo^s)=wmean(ixo^s,mom_c(idim))/rho(ixo^s)
2414  call twofl_get_csound_c_idim(wmean,x,ixi^l,ixo^l,idim,csoundr)
2415  if(present(cmin)) then
2416  cmax(ixo^s,1)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2417  cmin(ixo^s,1)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2418  if(h_correction) then
2419  {do ix^db=ixomin^db,ixomax^db\}
2420  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2421  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2422  {end do\}
2423  end if
2424  else
2425  cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
2426  end if
2427  !neutrals
2428 
2429  call get_rhon_tot(wmean,x,ixi^l,ixo^l,rho)
2430  tmp1(ixo^s)=wmean(ixo^s,mom_n(idim))/rho(ixo^s)
2431  call twofl_get_csound_n(wmean,x,ixi^l,ixo^l,csoundr)
2432  if(present(cmin)) then
2433  cmax(ixo^s,2)=max(abs(tmp1(ixo^s))+csoundr(ixo^s),zero)
2434  cmin(ixo^s,2)=min(abs(tmp1(ixo^s))-csoundr(ixo^s),zero)
2435  if(h_correction) then
2436  {do ix^db=ixomin^db,ixomax^db\}
2437  cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,2)),hspeed(ix^d,2))
2438  cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,2)),hspeed(ix^d,2))
2439  {end do\}
2440  end if
2441  else
2442  cmax(ixo^s,2)= abs(tmp1(ixo^s))+csoundr(ixo^s)
2443  end if
2444  case (3)
2445  ! Miyoshi 2005 JCP 208, 315 equation (67)
2446  call twofl_get_csound_c_idim(wlp,x,ixi^l,ixo^l,idim,csoundl)
2447  call twofl_get_csound_c_idim(wrp,x,ixi^l,ixo^l,idim,csoundr)
2448  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2449  if(present(cmin)) then
2450  cmin(ixo^s,1)=min(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))-csoundl(ixo^s)
2451  cmax(ixo^s,1)=max(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))+csoundl(ixo^s)
2452  if(h_correction) then
2453  {do ix^db=ixomin^db,ixomax^db\}
2454  cmin(ix^d,1)=sign(one,cmin(ix^d,1))*max(abs(cmin(ix^d,1)),hspeed(ix^d,1))
2455  cmax(ix^d,1)=sign(one,cmax(ix^d,1))*max(abs(cmax(ix^d,1)),hspeed(ix^d,1))
2456  {end do\}
2457  end if
2458  else
2459  cmax(ixo^s,1)=max(wlp(ixo^s,mom_c(idim)),wrp(ixo^s,mom_c(idim)))+csoundl(ixo^s)
2460  end if
2461  call twofl_get_csound_n(wlp,x,ixi^l,ixo^l,csoundl)
2462  call twofl_get_csound_n(wrp,x,ixi^l,ixo^l,csoundr)
2463  csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
2464  if(present(cmin)) then
2465  cmin(ixo^s,2)=min(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))-csoundl(ixo^s)
2466  cmax(ixo^s,2)=max(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))+csoundl(ixo^s)
2467  if(h_correction) then
2468  {do ix^db=ixomin^db,ixomax^db\}
2469  cmin(ix^d,2)=sign(one,cmin(ix^d,2))*max(abs(cmin(ix^d,1)),hspeed(ix^d,2))
2470  cmax(ix^d,2)=sign(one,cmax(ix^d,2))*max(abs(cmax(ix^d,1)),hspeed(ix^d,2))
2471  {end do\}
2472  end if
2473  else
2474  cmax(ixo^s,2)=max(wlp(ixo^s,mom_n(idim)),wrp(ixo^s,mom_n(idim)))+csoundl(ixo^s)
2475  end if
2476 
2477  end select
2478 
2479  end subroutine twofl_get_cbounds_species
2480 
2481  !> prepare velocities for ct methods
2482  subroutine twofl_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
2484 
2485  integer, intent(in) :: ixI^L, ixO^L, idim
2486  double precision, intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
2487  double precision, intent(in) :: cmax(ixI^S)
2488  double precision, intent(in), optional :: cmin(ixI^S)
2489  type(ct_velocity), intent(inout):: vcts
2490 
2491  integer :: idimE,idimN
2492 
2493  ! calculate velocities related to different UCT schemes
2494  select case(type_ct)
2495  case('average')
2496  case('uct_contact')
2497  if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
2498  ! get average normal velocity at cell faces
2499  vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom_c(idim))+wrp(ixo^s,mom_c(idim)))
2500  case('uct_hll')
2501  if(.not.allocated(vcts%vbarC)) then
2502  allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
2503  allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
2504  end if
2505  ! Store magnitude of characteristics
2506  if(present(cmin)) then
2507  vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
2508  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2509  else
2510  vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
2511  vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
2512  end if
2513 
2514  idimn=mod(idim,ndir)+1 ! 'Next' direction
2515  idime=mod(idim+1,ndir)+1 ! Electric field direction
2516  ! Store velocities
2517  vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom_c(idimn))
2518  vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom_c(idimn))
2519  vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
2520  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2521  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2522 
2523  vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom_c(idime))
2524  vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom_c(idime))
2525  vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
2526  +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
2527  /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
2528  case default
2529  call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
2530  end select
2531 
2532  end subroutine twofl_get_ct_velocity
2533 
2534  subroutine twofl_get_csound_c_idim(w,x,ixI^L,ixO^L,idim,csound)
2536 
2537  integer, intent(in) :: ixI^L, ixO^L, idim
2538  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2539  double precision, intent(out):: csound(ixI^S)
2540  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2541  double precision :: inv_rho(ixO^S)
2542  double precision :: tmp(ixI^S)
2543 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2544  double precision :: rhon(ixI^S)
2545 #endif
2546  call get_rhoc_tot(w,x,ixi^l,ixo^l,tmp)
2547 #if (!defined(ONE_FLUID) || ONE_FLUID==0) && (defined(A_TOT) && A_TOT == 1)
2548  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2549  inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+tmp(ixo^s))
2550 #else
2551  inv_rho(ixo^s)=1.d0/tmp(ixo^s)
2552 #endif
2553 
2554  call twofl_get_csound2_c_from_conserved(w,x,ixi^l,ixo^l,csound)
2555 
2556  ! store |B|^2 in v
2557  b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2558 
2559  cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2560  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2561  * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2562  * inv_rho(ixo^s)
2563 
2564  where(avmincs2(ixo^s)<zero)
2565  avmincs2(ixo^s)=zero
2566  end where
2567 
2568  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2569 
2570  if (.not. twofl_hall) then
2571  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2572  else
2573  ! take the Hall velocity into account:
2574  ! most simple estimate, high k limit:
2575  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2576  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2577  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2578  twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2579  end if
2580 
2581  end subroutine twofl_get_csound_c_idim
2582 
2583  !> Calculate fast magnetosonic wave speed when cbounds_species=false
2584  subroutine twofl_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
2586 
2587  integer, intent(in) :: ixI^L, ixO^L, idim
2588  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2589  double precision, intent(out):: csound(ixI^S)
2590  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2591  double precision :: inv_rho(ixO^S)
2592  double precision :: rhoc(ixI^S)
2593 #if (defined(A_TOT) && A_TOT == 1)
2594  double precision :: rhon(ixI^S)
2595 #endif
2596  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
2597 #if (defined(A_TOT) && A_TOT == 1)
2598  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2599  inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2600 #else
2601  inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2602 #endif
2603 
2604  call twofl_get_csound2_primitive(w,x,ixi^l,ixo^l,csound)
2605 
2606  ! store |B|^2 in v
2607  b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2608  cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2609  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2610  * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2611  * inv_rho(ixo^s)
2612 
2613  where(avmincs2(ixo^s)<zero)
2614  avmincs2(ixo^s)=zero
2615  end where
2616 
2617  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2618 
2619  if (.not. twofl_hall) then
2620  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2621  else
2622  ! take the Hall velocity into account:
2623  ! most simple estimate, high k limit:
2624  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2625  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2626  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2627  twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2628  end if
2629 
2630  contains
2631  !TODO copy it inside
2632  subroutine twofl_get_csound2_primitive(w,x,ixI^L,ixO^L,csound2)
2634  integer, intent(in) :: ixI^L, ixO^L
2635  double precision, intent(in) :: w(ixI^S,nw)
2636  double precision, intent(in) :: x(ixI^S,1:ndim)
2637  double precision, intent(out) :: csound2(ixI^S)
2638  double precision :: pth_c(ixI^S)
2639  double precision :: pth_n(ixI^S)
2640 
2641  if(phys_energy) then
2642  call twofl_get_pthermal_c_primitive(w,x,ixi^l,ixo^l,pth_c)
2643  call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,pth_n)
2644  call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2645  else
2646  call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2647  endif
2648  end subroutine twofl_get_csound2_primitive
2649 
2650  end subroutine twofl_get_csound_prim
2651 
2652  subroutine twofl_get_csound2(w,x,ixI^L,ixO^L,csound2)
2654  integer, intent(in) :: ixI^L, ixO^L
2655  double precision, intent(in) :: w(ixI^S,nw)
2656  double precision, intent(in) :: x(ixI^S,1:ndim)
2657  double precision, intent(out) :: csound2(ixI^S)
2658  double precision :: pth_c(ixI^S)
2659  double precision :: pth_n(ixI^S)
2660 
2661  if(phys_energy) then
2662  call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,pth_c)
2663  call twofl_get_pthermal_n(w,x,ixi^l,ixo^l,pth_n)
2664  call twofl_get_csound2_from_pthermal(w,x,ixi^l,ixo^l,pth_c,pth_n,csound2)
2665  else
2666  call twofl_get_csound2_adiab(w,x,ixi^l,ixo^l,csound2)
2667  endif
2668  end subroutine twofl_get_csound2
2669 
2670  subroutine twofl_get_csound2_adiab(w,x,ixI^L,ixO^L,csound2)
2672  integer, intent(in) :: ixI^L, ixO^L
2673  double precision, intent(in) :: w(ixI^S,nw)
2674  double precision, intent(in) :: x(ixI^S,1:ndim)
2675  double precision, intent(out) :: csound2(ixI^S)
2676  double precision :: rhoc(ixI^S)
2677  double precision :: rhon(ixI^S)
2678 
2679  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
2680  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2681  csound2(ixo^s)=twofl_gamma*twofl_adiab*&
2682  max((rhoc(ixo^s)**twofl_gamma + rhon(ixo^s)**twofl_gamma)/(rhoc(ixo^s)+ rhon(ixo^s)),&
2683  rhon(ixo^s)**gamma_1,rhoc(ixo^s)**gamma_1)
2684  end subroutine twofl_get_csound2_adiab
2685 
2686  subroutine twofl_get_csound(w,x,ixI^L,ixO^L,idim,csound)
2688 
2689  integer, intent(in) :: ixI^L, ixO^L, idim
2690  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2691  double precision, intent(out):: csound(ixI^S)
2692  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
2693  double precision :: inv_rho(ixO^S)
2694  double precision :: rhoc(ixI^S)
2695 #if (defined(A_TOT) && A_TOT == 1)
2696  double precision :: rhon(ixI^S)
2697 #endif
2698  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
2699 #if (defined(A_TOT) && A_TOT == 1)
2700  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2701  inv_rho(ixo^s) = 1d0/(rhon(ixo^s)+rhoc(ixo^s))
2702 #else
2703  inv_rho(ixo^s)=1.d0/rhoc(ixo^s)
2704 #endif
2705 
2706  call twofl_get_csound2(w,x,ixi^l,ixo^l,csound)
2707 
2708  ! store |B|^2 in v
2709  b2(ixo^s) = twofl_mag_en_all(w,ixi^l,ixo^l)
2710 
2711  cfast2(ixo^s) = b2(ixo^s) * inv_rho(ixo^s)+csound(ixo^s)
2712  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
2713  * twofl_mag_i_all(w,ixi^l,ixo^l,idim)**2 &
2714  * inv_rho(ixo^s)
2715 
2716  where(avmincs2(ixo^s)<zero)
2717  avmincs2(ixo^s)=zero
2718  end where
2719 
2720  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
2721 
2722  if (.not. twofl_hall) then
2723  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
2724  else
2725  ! take the Hall velocity into account:
2726  ! most simple estimate, high k limit:
2727  ! largest wavenumber supported by grid: Nyquist (in practise can reduce by some factor)
2728  kmax = dpi/min({dxlevel(^d)},bigdouble)*half
2729  csound(ixo^s) = max(sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s))), &
2730  twofl_etah * sqrt(b2(ixo^s))*inv_rho(ixo^s)*kmax)
2731  end if
2732 
2733  end subroutine twofl_get_csound
2734 
2735  subroutine twofl_get_csound2_from_pthermal(w,x,ixI^L,ixO^L,pth_c,pth_n,csound2)
2737  integer, intent(in) :: ixI^L, ixO^L
2738  double precision, intent(in) :: w(ixI^S,nw)
2739  double precision, intent(in) :: x(ixI^S,1:ndim)
2740  double precision, intent(in) :: pth_c(ixI^S)
2741  double precision, intent(in) :: pth_n(ixI^S)
2742  double precision, intent(out) :: csound2(ixI^S)
2743  double precision :: csound1(ixI^S),rhon(ixI^S),rhoc(ixI^S)
2744 
2745  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2746  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
2747 #if !defined(C_TOT) || C_TOT == 0
2748  csound2(ixo^s)=twofl_gamma*max((pth_c(ixo^s) + pth_n(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s)),&
2749  pth_n(ixo^s)/rhon(ixo^s), pth_c(ixo^s)/rhoc(ixo^s))
2750 #else
2751  csound2(ixo^s)=twofl_gamma*(csound2(ixo^s) + csound1(ixo^s))/(rhoc(ixo^s) + rhon(ixo^s))
2752 
2753 #endif
2754  end subroutine twofl_get_csound2_from_pthermal
2755 
2756 ! end cbounds_species=false
2757 
2758  subroutine twofl_get_csound_n(w,x,ixI^L,ixO^L,csound)
2760 
2761  integer, intent(in) :: ixI^L, ixO^L
2762  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
2763  double precision, intent(out):: csound(ixI^S)
2764  double precision :: pe_n1(ixI^S)
2765  call twofl_get_csound2_n_from_conserved(w,x,ixi^l,ixo^l,csound)
2766  csound(ixo^s) = sqrt(csound(ixo^s))
2767  end subroutine twofl_get_csound_n
2768 
2769  !> separate routines so that it is faster
2770  !> Calculate temperature=p/rho when in e_ the internal energy is stored
2771  subroutine twofl_get_temperature_from_eint_n(w, x, ixI^L, ixO^L, res)
2773  integer, intent(in) :: ixI^L, ixO^L
2774  double precision, intent(in) :: w(ixI^S, 1:nw)
2775  double precision, intent(in) :: x(ixI^S, 1:ndim)
2776  double precision, intent(out):: res(ixI^S)
2777 
2778  res(ixo^s) = 1d0/rn * gamma_1 * w(ixo^s, e_n_) /w(ixo^s,rho_n_)
2779 
2780  end subroutine twofl_get_temperature_from_eint_n
2781 
2782  subroutine twofl_get_temperature_from_eint_n_with_equi(w, x, ixI^L, ixO^L, res)
2784  integer, intent(in) :: ixI^L, ixO^L
2785  double precision, intent(in) :: w(ixI^S, 1:nw)
2786  double precision, intent(in) :: x(ixI^S, 1:ndim)
2787  double precision, intent(out):: res(ixI^S)
2788 
2789  res(ixo^s) = 1d0/rn * (gamma_1 * w(ixo^s, e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)) /&
2790  (w(ixo^s,rho_n_) +block%equi_vars(ixo^s,equi_rho_n0_,b0i))
2792 
2793 ! subroutine twofl_get_temperature_n_pert_from_tot(Te, ixI^L, ixO^L, res)
2794 ! use mod_global_parameters
2795 ! integer, intent(in) :: ixI^L, ixO^L
2796 ! double precision, intent(in) :: Te(ixI^S)
2797 ! double precision, intent(out):: res(ixI^S)
2798 ! res(ixO^S) = Te(ixO^S) -1d0/Rn * &
2799 ! block%equi_vars(ixO^S,equi_pe_n0_,0)/block%equi_vars(ixO^S,equi_rho_n0_,0)
2800 ! end subroutine twofl_get_temperature_n_pert_from_tot
2801 
2802  subroutine twofl_get_temperature_n_equi(w,x, ixI^L, ixO^L, res)
2804  integer, intent(in) :: ixI^L, ixO^L
2805  double precision, intent(in) :: w(ixI^S, 1:nw)
2806  double precision, intent(in) :: x(ixI^S, 1:ndim)
2807  double precision, intent(out):: res(ixI^S)
2808  res(ixo^s) = 1d0/rn * &
2809  block%equi_vars(ixo^s,equi_pe_n0_,b0i)/block%equi_vars(ixo^s,equi_rho_n0_,b0i)
2810  end subroutine twofl_get_temperature_n_equi
2811 
2812  subroutine twofl_get_rho_n_equi(w, x,ixI^L, ixO^L, res)
2814  integer, intent(in) :: ixI^L, ixO^L
2815  double precision, intent(in) :: w(ixI^S, 1:nw)
2816  double precision, intent(in) :: x(ixI^S, 1:ndim)
2817  double precision, intent(out):: res(ixI^S)
2818  res(ixo^s) = block%equi_vars(ixo^s,equi_rho_n0_,b0i)
2819  end subroutine twofl_get_rho_n_equi
2820 
2821  subroutine twofl_get_pe_n_equi(w, x, ixI^L, ixO^L, res)
2823  integer, intent(in) :: ixI^L, ixO^L
2824  double precision, intent(in) :: w(ixI^S, 1:nw)
2825  double precision, intent(in) :: x(ixI^S, 1:ndim)
2826  double precision, intent(out):: res(ixI^S)
2827  res(ixo^s) = block%equi_vars(ixo^s,equi_pe_n0_,b0i)
2828  end subroutine twofl_get_pe_n_equi
2829 
2830  !> Calculate temperature=p/rho when in e_ the total energy is stored
2831  !> this does not check the values of twofl_energy and twofl_internal_e,
2832  !> twofl_energy = .true. and twofl_internal_e = .false.
2833  !> also check small_values is avoided
2834  subroutine twofl_get_temperature_from_etot_n(w, x, ixI^L, ixO^L, res)
2836  integer, intent(in) :: ixI^L, ixO^L
2837  double precision, intent(in) :: w(ixI^S, 1:nw)
2838  double precision, intent(in) :: x(ixI^S, 1:ndim)
2839  double precision, intent(out):: res(ixI^S)
2840  res(ixo^s)=1d0/rn * (gamma_1*(w(ixo^s,e_n_)&
2841  - twofl_kin_en_n(w,ixi^l,ixo^l)))/w(ixo^s,rho_n_)
2842  end subroutine twofl_get_temperature_from_etot_n
2843 
2844  subroutine twofl_get_temperature_from_etot_n_with_equi(w, x, ixI^L, ixO^L, res)
2846  integer, intent(in) :: ixI^L, ixO^L
2847  double precision, intent(in) :: w(ixI^S, 1:nw)
2848  double precision, intent(in) :: x(ixI^S, 1:ndim)
2849  double precision, intent(out):: res(ixI^S)
2850  res(ixo^s)=1d0/rn * (gamma_1*(w(ixo^s,e_n_)&
2851  - twofl_kin_en_n(w,ixi^l,ixo^l)) + block%equi_vars(ixo^s,equi_pe_n0_,b0i))&
2852  /(w(ixo^s,rho_n_) +block%equi_vars(ixo^s,equi_rho_n0_,b0i))
2853 
2855 
2856  !> separate routines so that it is faster
2857  !> Calculate temperature=p/rho when in e_ the internal energy is stored
2858  subroutine twofl_get_temperature_from_eint_c(w, x, ixI^L, ixO^L, res)
2860  integer, intent(in) :: ixI^L, ixO^L
2861  double precision, intent(in) :: w(ixI^S, 1:nw)
2862  double precision, intent(in) :: x(ixI^S, 1:ndim)
2863  double precision, intent(out):: res(ixI^S)
2864 
2865  res(ixo^s) = 1d0/rc * gamma_1 * w(ixo^s, e_c_) /w(ixo^s,rho_c_)
2866 
2867  end subroutine twofl_get_temperature_from_eint_c
2868 
2869  subroutine twofl_get_temperature_from_eint_c_with_equi(w, x, ixI^L, ixO^L, res)
2871  integer, intent(in) :: ixI^L, ixO^L
2872  double precision, intent(in) :: w(ixI^S, 1:nw)
2873  double precision, intent(in) :: x(ixI^S, 1:ndim)
2874  double precision, intent(out):: res(ixI^S)
2875  res(ixo^s) = 1d0/rc * (gamma_1 * w(ixo^s, e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)) /&
2876  (w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
2878 
2879 ! subroutine twofl_get_temperature_c_pert_from_tot(Te, ixI^L, ixO^L, res)
2880 ! use mod_global_parameters
2881 ! integer, intent(in) :: ixI^L, ixO^L
2882 ! double precision, intent(in) :: Te(ixI^S)
2883 ! double precision, intent(out):: res(ixI^S)
2884 ! res(ixO^S) = Te(ixO^S) -1d0/Rc * &
2885 ! block%equi_vars(ixO^S,equi_pe_c0_,0)/block%equi_vars(ixO^S,equi_rho_c0_,0)
2886 ! end subroutine twofl_get_temperature_c_pert_from_tot
2887 
2888  subroutine twofl_get_temperature_c_equi(w,x, ixI^L, ixO^L, res)
2890  integer, intent(in) :: ixI^L, ixO^L
2891  double precision, intent(in) :: w(ixI^S, 1:nw)
2892  double precision, intent(in) :: x(ixI^S, 1:ndim)
2893  double precision, intent(out):: res(ixI^S)
2894  res(ixo^s) = 1d0/rc * &
2895  block%equi_vars(ixo^s,equi_pe_c0_,b0i)/block%equi_vars(ixo^s,equi_rho_c0_,b0i)
2896  end subroutine twofl_get_temperature_c_equi
2897 
2898  subroutine twofl_get_rho_c_equi(w, x, ixI^L, ixO^L, res)
2900  integer, intent(in) :: ixI^L, ixO^L
2901  double precision, intent(in) :: w(ixI^S, 1:nw)
2902  double precision, intent(in) :: x(ixI^S, 1:ndim)
2903  double precision, intent(out):: res(ixI^S)
2904  res(ixo^s) = block%equi_vars(ixo^s,equi_rho_c0_,b0i)
2905  end subroutine twofl_get_rho_c_equi
2906 
2907  subroutine twofl_get_pe_c_equi(w,x, ixI^L, ixO^L, res)
2909  integer, intent(in) :: ixI^L, ixO^L
2910  double precision, intent(in) :: w(ixI^S, 1:nw)
2911  double precision, intent(in) :: x(ixI^S, 1:ndim)
2912  double precision, intent(out):: res(ixI^S)
2913  res(ixo^s) = block%equi_vars(ixo^s,equi_pe_c0_,b0i)
2914  end subroutine twofl_get_pe_c_equi
2915 
2916  !> Calculate temperature=p/rho when in e_ the total energy is stored
2917  !> this does not check the values of twofl_energy and twofl_internal_e,
2918  !> twofl_energy = .true. and twofl_internal_e = .false.
2919  !> also check small_values is avoided
2920  subroutine twofl_get_temperature_from_etot_c(w, x, ixI^L, ixO^L, res)
2922  integer, intent(in) :: ixI^L, ixO^L
2923  double precision, intent(in) :: w(ixI^S, 1:nw)
2924  double precision, intent(in) :: x(ixI^S, 1:ndim)
2925  double precision, intent(out):: res(ixI^S)
2926  res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
2927  - twofl_kin_en_c(w,ixi^l,ixo^l)&
2928  - twofl_mag_en(w,ixi^l,ixo^l)))/w(ixo^s,rho_c_)
2929  end subroutine twofl_get_temperature_from_etot_c
2930  subroutine twofl_get_temperature_from_eki_c(w, x, ixI^L, ixO^L, res)
2932  integer, intent(in) :: ixI^L, ixO^L
2933  double precision, intent(in) :: w(ixI^S, 1:nw)
2934  double precision, intent(in) :: x(ixI^S, 1:ndim)
2935  double precision, intent(out):: res(ixI^S)
2936  res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
2937  - twofl_kin_en_c(w,ixi^l,ixo^l)))/w(ixo^s,rho_c_)
2938  end subroutine twofl_get_temperature_from_eki_c
2939 
2940  subroutine twofl_get_temperature_from_etot_c_with_equi(w, x, ixI^L, ixO^L, res)
2942  integer, intent(in) :: ixI^L, ixO^L
2943  double precision, intent(in) :: w(ixI^S, 1:nw)
2944  double precision, intent(in) :: x(ixI^S, 1:ndim)
2945  double precision, intent(out):: res(ixI^S)
2946  res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
2947  - twofl_kin_en_c(w,ixi^l,ixo^l)&
2948  - twofl_mag_en(w,ixi^l,ixo^l)) + block%equi_vars(ixo^s,equi_pe_c0_,b0i))&
2949  /(w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
2950 
2952 
2953  subroutine twofl_get_temperature_from_eki_c_with_equi(w, x, ixI^L, ixO^L, res)
2955  integer, intent(in) :: ixI^L, ixO^L
2956  double precision, intent(in) :: w(ixI^S, 1:nw)
2957  double precision, intent(in) :: x(ixI^S, 1:ndim)
2958  double precision, intent(out):: res(ixI^S)
2959  res(ixo^s)=1d0/rc * (gamma_1*(w(ixo^s,e_c_)&
2960  - twofl_kin_en_c(w,ixi^l,ixo^l)) + block%equi_vars(ixo^s,equi_pe_c0_,b0i))&
2961  /(w(ixo^s,rho_c_) +block%equi_vars(ixo^s,equi_rho_c0_,b0i))
2962 
2964 
2965  subroutine twofl_get_csound2_adiab_n(w,x,ixI^L,ixO^L,csound2)
2967  integer, intent(in) :: ixI^L, ixO^L
2968  double precision, intent(in) :: w(ixI^S,nw)
2969  double precision, intent(in) :: x(ixI^S,1:ndim)
2970  double precision, intent(out) :: csound2(ixI^S)
2971  double precision :: rhon(ixI^S)
2972 
2973  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2974  csound2(ixo^s)=twofl_gamma*twofl_adiab*rhon(ixo^s)**gamma_1
2975 
2976  end subroutine twofl_get_csound2_adiab_n
2977 
2978  subroutine twofl_get_csound2_n_from_conserved(w,x,ixI^L,ixO^L,csound2)
2980  integer, intent(in) :: ixI^L, ixO^L
2981  double precision, intent(in) :: w(ixI^S,nw)
2982  double precision, intent(in) :: x(ixI^S,1:ndim)
2983  double precision, intent(out) :: csound2(ixI^S)
2984  double precision :: rhon(ixI^S)
2985 
2986  if(phys_energy) then
2987  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
2988  call twofl_get_pthermal_n(w,x,ixi^l,ixo^l,csound2)
2989  csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
2990  else
2991  call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
2992  endif
2993  end subroutine twofl_get_csound2_n_from_conserved
2994 
2995  !! TO DELETE
2996  subroutine twofl_get_csound2_n_from_primitive(w,x,ixI^L,ixO^L,csound2)
2998  integer, intent(in) :: ixI^L, ixO^L
2999  double precision, intent(in) :: w(ixI^S,nw)
3000  double precision, intent(in) :: x(ixI^S,1:ndim)
3001  double precision, intent(out) :: csound2(ixI^S)
3002  double precision :: rhon(ixI^S)
3003 
3004  if(phys_energy) then
3005  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
3006  call twofl_get_pthermal_n_primitive(w,x,ixi^l,ixo^l,csound2)
3007  csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhon(ixo^s)
3008  else
3009  call twofl_get_csound2_adiab_n(w,x,ixi^l,ixo^l,csound2)
3010  endif
3011  end subroutine twofl_get_csound2_n_from_primitive
3012 
3013  subroutine twofl_get_csound2_adiab_c(w,x,ixI^L,ixO^L,csound2)
3015  integer, intent(in) :: ixI^L, ixO^L
3016  double precision, intent(in) :: w(ixI^S,nw)
3017  double precision, intent(in) :: x(ixI^S,1:ndim)
3018  double precision, intent(out) :: csound2(ixI^S)
3019  double precision :: rhoc(ixI^S)
3020 
3021  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3022  csound2(ixo^s)=twofl_gamma*twofl_adiab* rhoc(ixo^s)**gamma_1
3023 
3024  end subroutine twofl_get_csound2_adiab_c
3025 
3026  subroutine twofl_get_csound2_c_from_conserved(w,x,ixI^L,ixO^L,csound2)
3028  integer, intent(in) :: ixi^l, ixo^l
3029  double precision, intent(in) :: w(ixi^s,nw)
3030  double precision, intent(in) :: x(ixi^s,1:ndim)
3031  double precision, intent(out) :: csound2(ixi^s)
3032  double precision :: rhoc(ixi^s)
3033 
3034  if(phys_energy) then
3035  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3036  call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,csound2)
3037  csound2(ixo^s)=twofl_gamma*csound2(ixo^s)/rhoc(ixo^s)
3038  else
3039  call twofl_get_csound2_adiab_c(w,x,ixi^l,ixo^l,csound2)
3040  endif
3041  end subroutine twofl_get_csound2_c_from_conserved
3042 
3043  !> Calculate fluxes within ixO^L.
3044  subroutine twofl_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
3046  use mod_geometry
3047 
3048  integer, intent(in) :: ixI^L, ixO^L, idim
3049  ! conservative w
3050  double precision, intent(in) :: wC(ixI^S,nw)
3051  ! primitive w
3052  double precision, intent(in) :: w(ixI^S,nw)
3053  double precision, intent(in) :: x(ixI^S,1:ndim)
3054  double precision,intent(out) :: f(ixI^S,nwflux)
3055 
3056  double precision :: pgas(ixO^S), ptotal(ixO^S),tmp(ixI^S)
3057  double precision, allocatable:: vHall(:^D&,:)
3058  integer :: idirmin, iw, idir, jdir, kdir
3059 
3060  ! value at the interfaces, idim = block%iw0 --> b0i
3061  ! reuse tmp, used afterwards
3062  ! value at the interface so we can't put momentum
3063  call get_rhoc_tot(w,x,ixi^l,ixo^l,tmp)
3064  ! Get flux of density
3065  f(ixo^s,rho_c_)=w(ixo^s,mom_c(idim))*tmp(ixo^s)
3066  ! pgas is time dependent only
3067  if(phys_energy) then
3068  pgas(ixo^s)=w(ixo^s,e_c_)
3069  else
3070  pgas(ixo^s)=twofl_adiab*tmp(ixo^s)**twofl_gamma
3071  if(has_equi_pe_c0) then
3072  pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe_c0_,b0i)
3073  end if
3074  end if
3075 
3076  if (twofl_hall) then
3077  allocate(vhall(ixi^s,1:ndir))
3078  call twofl_getv_hall(w,x,ixi^l,ixo^l,vhall)
3079  end if
3080 
3081  if(b0field) tmp(ixo^s)=sum(block%B0(ixo^s,:,idim)*w(ixo^s,mag(:)),dim=ndim+1)
3082 
3083  ptotal(ixo^s) = pgas(ixo^s) + 0.5d0*sum(w(ixo^s, mag(:))**2, dim=ndim+1)
3084 
3085  ! Get flux of momentum
3086  ! f_i[m_k]=v_i*m_k-b_k*b_i [+ptotal if i==k]
3087  do idir=1,ndir
3088  if(idim==idir) then
3089  f(ixo^s,mom_c(idir))=ptotal(ixo^s)-w(ixo^s,mag(idim))*w(ixo^s,mag(idir))
3090  if(b0field) f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))+tmp(ixo^s)
3091  else
3092  f(ixo^s,mom_c(idir))= -w(ixo^s,mag(idir))*w(ixo^s,mag(idim))
3093  end if
3094  if (b0field) then
3095  f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))&
3096  -w(ixo^s,mag(idir))*block%B0(ixo^s,idim,idim)&
3097  -w(ixo^s,mag(idim))*block%B0(ixo^s,idir,idim)
3098  end if
3099  f(ixo^s,mom_c(idir))=f(ixo^s,mom_c(idir))+w(ixo^s,mom_c(idim))*wc(ixo^s,mom_c(idir))
3100  end do
3101 
3102  ! Get flux of energy
3103  ! f_i[e]=v_i*e+v_i*ptotal-b_i*(b_k*v_k)
3104  if(phys_energy) then
3105  if (phys_internal_e) then
3106  f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*wc(ixo^s,e_c_)
3107  else if(twofl_eq_energy == eq_energy_ki) then
3108 
3109  f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*(wc(ixo^s,e_c_)+pgas(ixo^s))
3110  else
3111  f(ixo^s,e_c_)=w(ixo^s,mom_c(idim))*(wc(ixo^s,e_c_)+ptotal(ixo^s))&
3112  -w(ixo^s,mag(idim))*sum(w(ixo^s,mag(:))*w(ixo^s,mom_c(:)),dim=ndim+1)
3113 
3114  if (b0field) then
3115  f(ixo^s,e_c_) = f(ixo^s,e_c_) &
3116  + w(ixo^s,mom_c(idim)) * tmp(ixo^s) &
3117  - sum(w(ixo^s,mom_c(:))*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
3118  end if
3119 
3120  if (twofl_hall) then
3121  ! f_i[e]= f_i[e] + vHall_i*(b_k*b_k) - b_i*(vHall_k*b_k)
3122  if (twofl_etah>zero) then
3123  f(ixo^s,e_c_) = f(ixo^s,e_c_) + vhall(ixo^s,idim) * &
3124  sum(w(ixo^s, mag(:))**2,dim=ndim+1) &
3125  - w(ixo^s,mag(idim)) * sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1)
3126  if (b0field) then
3127  f(ixo^s,e_c_) = f(ixo^s,e_c_) &
3128  + vhall(ixo^s,idim) * tmp(ixo^s) &
3129  - sum(vhall(ixo^s,:)*w(ixo^s,mag(:)),dim=ndim+1) * block%B0(ixo^s,idim,idim)
3130  end if
3131  end if
3132  end if
3133  end if !total_energy
3134  ! add flux of equilibrium internal energy corresponding to pe_c0
3135  if(has_equi_pe_c0) then
3136 #if !defined(E_RM_W0) || E_RM_W0 == 1
3137  f(ixo^s,e_c_)= f(ixo^s,e_c_) &
3138  + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * inv_gamma_1
3139 #else
3140  if(phys_internal_e) then
3141  f(ixo^s,e_c_)= f(ixo^s,e_c_) &
3142  + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * inv_gamma_1
3143  else
3144  f(ixo^s,e_c_)= f(ixo^s,e_c_) &
3145  + w(ixo^s,mom_c(idim)) * block%equi_vars(ixo^s,equi_pe_c0_,idim) * twofl_gamma * inv_gamma_1
3146  end if
3147 #endif
3148  end if
3149  end if !phys_energy
3150 
3151  ! compute flux of magnetic field
3152  ! f_i[b_k]=v_i*b_k-v_k*b_i
3153  do idir=1,ndir
3154  if (idim==idir) then
3155  ! f_i[b_i] should be exactly 0, so we do not use the transport flux
3156  if (twofl_glm) then
3157  f(ixo^s,mag(idir))=w(ixo^s,psi_)
3158  else
3159  f(ixo^s,mag(idir))=zero
3160  end if
3161  else
3162  f(ixo^s,mag(idir))=w(ixo^s,mom_c(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom_c(idir))
3163 
3164  if (b0field) then
3165  f(ixo^s,mag(idir))=f(ixo^s,mag(idir))&
3166  +w(ixo^s,mom_c(idim))*block%B0(ixo^s,idir,idim)&
3167  -w(ixo^s,mom_c(idir))*block%B0(ixo^s,idim,idim)
3168  end if
3169 
3170  if (twofl_hall) then
3171  ! f_i[b_k] = f_i[b_k] + vHall_i*b_k - vHall_k*b_i
3172  if (twofl_etah>zero) then
3173  if (b0field) then
3174  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3175  - vhall(ixo^s,idir)*(w(ixo^s,mag(idim))+block%B0(ixo^s,idim,idim)) &
3176  + vhall(ixo^s,idim)*(w(ixo^s,mag(idir))+block%B0(ixo^s,idir,idim))
3177  else
3178  f(ixo^s,mag(idir)) = f(ixo^s,mag(idir)) &
3179  - vhall(ixo^s,idir)*w(ixo^s,mag(idim)) &
3180  + vhall(ixo^s,idim)*w(ixo^s,mag(idir))
3181  end if
3182  end if
3183  end if
3184 
3185  end if
3186  end do
3187 
3188  if (twofl_glm) then
3189  !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
3190  f(ixo^s,psi_) = cmax_global**2*w(ixo^s,mag(idim))
3191  end if
3192 
3193  if (twofl_hall) then
3194  deallocate(vhall)
3195  end if
3196 
3197  !!neutrals
3198  call get_rhon_tot(w,x,ixi^l,ixo^l,tmp)
3199  f(ixo^s,rho_n_)=w(ixo^s,mom_n(idim))*tmp(ixo^s)
3200  if(phys_energy) then
3201  pgas(ixo^s) = w(ixo^s, e_n_)
3202  else
3203  pgas(ixo^s)=twofl_adiab*tmp(ixo^s)**twofl_gamma
3204  if(has_equi_pe_n0) then
3205  pgas(ixo^s)=pgas(ixo^s)-block%equi_vars(ixo^s,equi_pe_n0_,b0i)
3206  end if
3207  end if
3208  ! Momentum flux is v_i*m_i, +p in direction idim
3209  do idir = 1, ndir
3210  !if(idim==idir) then
3211  ! f(ixO^S,mom_c(idir)) = pgas(ixO^S)
3212  !else
3213  ! f(ixO^S,mom_c(idir)) = 0.0d0
3214  !end if
3215  !f(ixO^S,mom_c(idir))=f(ixO^S,mom_c(idir))+w(ixO^S,mom_c(idim))*wC(ixO^S,mom_c(idir))
3216  f(ixo^s, mom_n(idir)) = w(ixo^s,mom_n(idim)) * wc(ixo^s, mom_n(idir))
3217  end do
3218 
3219  f(ixo^s, mom_n(idim)) = f(ixo^s, mom_n(idim)) + pgas(ixo^s)
3220 
3221  if(phys_energy) then
3222  !reuse pgas for storing a in the term: div (u_n * a) and make multiplication at the end
3223  pgas(ixo^s) = wc(ixo^s,e_n_)
3224  if(.not. phys_internal_e) then
3225  ! add pressure perturbation
3226  pgas(ixo^s) = pgas(ixo^s) + w(ixo^s,e_n_)
3227  end if
3228  ! add flux of equilibrium internal energy corresponding to pe_n0
3229  if(has_equi_pe_n0) then
3230 #if !defined(E_RM_W0) || E_RM_W0 == 1
3231  pgas(ixo^s) = pgas(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,idim) * inv_gamma_1
3232 #else
3233  pgas(ixo^s) = pgas(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,idim) * twofl_gamma * inv_gamma_1
3234 #endif
3235  end if
3236  ! add u_n * a in the flux
3237  f(ixo^s, e_n_) = w(ixo^s,mom_n(idim)) * pgas(ixo^s)
3238 
3239  ! Viscosity fluxes - viscInDiv
3240  !if (hd_viscosity) then
3241  ! call visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, phys_energy)
3242  !endif
3243  end if
3244 
3245  end subroutine twofl_get_flux
3246 
3247  !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
3248  subroutine twofl_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
3252  !use mod_gravity, only: gravity_add_source
3253 
3254  integer, intent(in) :: ixI^L, ixO^L
3255  double precision, intent(in) :: qdt,dtfactor
3256  double precision, intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw),x(ixI^S,1:ndim)
3257  double precision, intent(inout) :: w(ixI^S,1:nw)
3258  logical, intent(in) :: qsourcesplit
3259  logical, intent(inout) :: active
3260 
3261  if (.not. qsourcesplit) then
3262  ! Source for solving internal energy
3263  if(phys_internal_e) then
3264  active = .true.
3265  call internal_energy_add_source_n(qdt,ixi^l,ixo^l,wct,w,x)
3266  call internal_energy_add_source_c(qdt,ixi^l,ixo^l,wct,w,x,e_c_)
3267  else
3268 #if !defined(E_RM_W0) || E_RM_W0==1
3269  ! add -p0 div v source terms when equi are present
3270  if(has_equi_pe_n0) then
3271  active = .true.
3272  call add_pe_n0_divv(qdt,ixi^l,ixo^l,wct,w,x)
3273  endif
3274  if(has_equi_pe_c0) then
3275  active = .true.
3276  call add_pe_c0_divv(qdt,ixi^l,ixo^l,wct,w,x)
3277  endif
3278 #endif
3279  if(twofl_eq_energy == eq_energy_ki) then
3280  active = .true.
3281  call add_source_lorentz_work(qdt,ixi^l,ixo^l,w,wct,x)
3282  endif
3283  endif
3284 
3285  ! Source for B0 splitting
3286  if (b0field) then
3287  active = .true.
3288  call add_source_b0split(qdt,ixi^l,ixo^l,wct,w,x)
3289  end if
3290 
3291  ! Sources for resistivity in eqs. for e, B1, B2 and B3
3292  if (abs(twofl_eta)>smalldouble)then
3293  active = .true.
3294  call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
3295  end if
3296 
3297  if (twofl_eta_hyper>0.d0)then
3298  active = .true.
3299  call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
3300  end if
3301  !it is not added in a split manner
3302  if(.not. use_imex_scheme .and. has_collisions()) then
3303  active = .true.
3304  call twofl_explicit_coll_terms_update(qdt,ixi^l,ixo^l,w,wct,x)
3305  endif
3306 
3307  if(twofl_hyperdiffusivity) then
3308  active = .true.
3309  call add_source_hyperdiffusive(qdt,ixi^l,ixo^l,w,wct,x)
3310  endif
3311 
3312  end if
3313 
3314  {^nooned
3315  if(source_split_divb .eqv. qsourcesplit) then
3316  ! Sources related to div B
3317  select case (type_divb)
3318  case (divb_none)
3319  ! Do nothing
3320  case (divb_glm)
3321  active = .true.
3322  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
3323  case (divb_powel)
3324  active = .true.
3325  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
3326  case (divb_janhunen)
3327  active = .true.
3328  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
3329  case (divb_linde)
3330  active = .true.
3331  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
3332  case (divb_lindejanhunen)
3333  active = .true.
3334  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
3335  call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
3336  case (divb_lindepowel)
3337  active = .true.
3338  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
3339  call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
3340  case (divb_lindeglm)
3341  active = .true.
3342  call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
3343  call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
3344  case (divb_ct)
3345  continue ! Do nothing
3346  case (divb_multigrid)
3347  continue ! Do nothing
3348  case default
3349  call mpistop('Unknown divB fix')
3350  end select
3351  end if
3352  }
3353 
3354  if(twofl_radiative_cooling_c) then
3355  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
3356  w,x,qsourcesplit,active,rc_fl_c)
3357  end if
3358  if(twofl_radiative_cooling_n) then
3359  call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
3360  w,x,qsourcesplit,active,rc_fl_n)
3361  end if
3362 !
3363 ! if(twofl_viscosity) then
3364 ! call viscosity_add_source(qdt,ixI^L,ixO^L,wCT,&
3365 ! w,x,phys_energy,qsourcesplit,active)
3366 ! end if
3367 !
3368  if(twofl_gravity) then
3369  call gravity_add_source(qdt,ixi^l,ixo^l,wct,&
3370  w,x,twofl_eq_energy .eq. eq_energy_ki .or. phys_total_energy,qsourcesplit,active)
3371  end if
3372 
3373  end subroutine twofl_add_source
3374 
3375  subroutine add_pe_n0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3377  use mod_geometry
3378 
3379  integer, intent(in) :: ixI^L, ixO^L
3380  double precision, intent(in) :: qdt
3381  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3382  double precision, intent(inout) :: w(ixI^S,1:nw)
3383  double precision :: v(ixI^S,1:ndir)
3384 
3385  call twofl_get_v_n(wct,x,ixi^l,ixi^l,v)
3386  call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-block%equi_vars(ixi^s,equi_pe_n0_,0),w,x,e_n_)
3387 
3388  end subroutine add_pe_n0_divv
3389 
3390  subroutine add_pe_c0_divv(qdt,ixI^L,ixO^L,wCT,w,x)
3392  use mod_geometry
3393 
3394  integer, intent(in) :: ixI^L, ixO^L
3395  double precision, intent(in) :: qdt
3396  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3397  double precision, intent(inout) :: w(ixI^S,1:nw)
3398  double precision :: v(ixI^S,1:ndir)
3399 
3400  call twofl_get_v_c(wct,x,ixi^l,ixi^l,v)
3401  call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-block%equi_vars(ixi^s,equi_pe_c0_,0),w,x,e_c_)
3402 
3403  end subroutine add_pe_c0_divv
3404 
3405  subroutine add_geom_pdivv(qdt,ixI^L,ixO^L,v,p,w,x,ind)
3407  use mod_geometry
3408 
3409  integer, intent(in) :: ixI^L, ixO^L,ind
3410  double precision, intent(in) :: qdt
3411  double precision, intent(in) :: p(ixI^S), v(ixI^S,1:ndir), x(ixI^S,1:ndim)
3412  double precision, intent(inout) :: w(ixI^S,1:nw)
3413  double precision :: divv(ixI^S)
3414 
3415  if(slab_uniform) then
3416  if(nghostcells .gt. 2) then
3417  call divvector(v,ixi^l,ixo^l,divv,sixthorder=.true.)
3418  else
3419  call divvector(v,ixi^l,ixo^l,divv,fourthorder=.true.)
3420  end if
3421  else
3422  call divvector(v,ixi^l,ixo^l,divv)
3423  end if
3424  w(ixo^s,ind)=w(ixo^s,ind)+qdt*p(ixo^s)*divv(ixo^s)
3425  end subroutine add_geom_pdivv
3426 
3427  !> Compute the Lorentz force (JxB)
3428  subroutine get_lorentz(ixI^L,ixO^L,w,JxB)
3430  integer, intent(in) :: ixI^L, ixO^L
3431  double precision, intent(in) :: w(ixI^S,1:nw)
3432  double precision, intent(inout) :: JxB(ixI^S,3)
3433  double precision :: a(ixI^S,3), b(ixI^S,3), tmp(ixI^S,3)
3434  integer :: idir, idirmin
3435  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
3436  double precision :: current(ixI^S,7-2*ndir:3)
3437 
3438  b=0.0d0
3439  do idir = 1, ndir
3440  b(ixo^s, idir) = twofl_mag_i_all(w, ixi^l, ixo^l,idir)
3441  end do
3442 
3443  ! store J current in a
3444  call get_current(w,ixi^l,ixo^l,idirmin,current)
3445 
3446  a=0.0d0
3447  do idir=7-2*ndir,3
3448  a(ixo^s,idir)=current(ixo^s,idir)
3449  end do
3450 
3451  call cross_product(ixi^l,ixo^l,a,b,jxb)
3452  end subroutine get_lorentz
3453 
3454  subroutine add_source_lorentz_work(qdt,ixI^L,ixO^L,w,wCT,x)
3456  integer, intent(in) :: ixI^L, ixO^L
3457  double precision, intent(in) :: qdt
3458  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3459  double precision, intent(inout) :: w(ixI^S,1:nw)
3460  double precision :: a(ixI^S,3), b(ixI^S,1:ndir)
3461 
3462  call get_lorentz(ixi^l, ixo^l,wct,a)
3463  call twofl_get_v_c(wct,x,ixi^l,ixo^l,b)
3464  w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*sum(a(ixo^s,1:ndir)*b(ixo^s,1:ndir),dim=ndim+1)
3465 
3466  end subroutine add_source_lorentz_work
3467 
3468  !> Calculate v_n vector
3469  subroutine twofl_get_v_n(w,x,ixI^L,ixO^L,v)
3471 
3472  integer, intent(in) :: ixI^L, ixO^L
3473  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3474  double precision, intent(out) :: v(ixI^S,ndir)
3475  double precision :: rhon(ixI^S)
3476  integer :: idir
3477 
3478  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
3479 
3480  do idir=1,ndir
3481  v(ixo^s,idir) = w(ixo^s, mom_n(idir)) / rhon(ixo^s)
3482  end do
3483 
3484  end subroutine twofl_get_v_n
3485 
3486  subroutine get_rhon_tot(w,x,ixI^L,ixO^L,rhon)
3488  integer, intent(in) :: ixi^l, ixo^l
3489  double precision, intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
3490  double precision, intent(out) :: rhon(ixi^s)
3491  if(has_equi_rho_n0) then
3492  rhon(ixo^s) = w(ixo^s,rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,b0i)
3493  else
3494  rhon(ixo^s) = w(ixo^s,rho_n_)
3495  endif
3496 
3497  end subroutine get_rhon_tot
3498 
3499  subroutine twofl_get_pthermal_n(w,x,ixI^L,ixO^L,pth)
3502  integer, intent(in) :: ixi^l, ixo^l
3503  double precision, intent(in) :: w(ixi^s,1:nw)
3504  double precision, intent(in) :: x(ixi^s,1:ndim)
3505  double precision, intent(out) :: pth(ixi^s)
3506 
3507  integer :: ix^d, iw
3508 
3509  if(phys_energy) then
3510  if(phys_internal_e) then
3511  pth(ixo^s)=gamma_1*w(ixo^s,e_n_)
3512  else
3513  pth(ixo^s)=gamma_1*(w(ixo^s,e_n_)&
3514  - twofl_kin_en_n(w,ixi^l,ixo^l))
3515  end if
3516  if(has_equi_pe_n0) then
3517  pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
3518  endif
3519  else
3520  call get_rhon_tot(w,x,ixi^l,ixo^l,pth)
3521  pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
3522  end if
3523 
3524  if (fix_small_values) then
3525  {do ix^db= ixo^lim^db\}
3526  if(pth(ix^d)<small_pressure) then
3527  pth(ix^d)=small_pressure
3528  end if
3529  {enddo^d&\}
3530  else if (check_small_values) then
3531  {do ix^db= ixo^lim^db\}
3532  if(pth(ix^d)<small_pressure) then
3533  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3534  " encountered when call twofl_get_pthermal_n"
3535  write(*,*) "Iteration: ", it, " Time: ", global_time
3536  write(*,*) "Location: ", x(ix^d,:)
3537  write(*,*) "Cell number: ", ix^d
3538  do iw=1,nw
3539  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3540  end do
3541  ! use erroneous arithmetic operation to crash the run
3542  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3543  write(*,*) "Saving status at the previous time step"
3544  crash=.true.
3545  end if
3546  {enddo^d&\}
3547  end if
3548 
3549  end subroutine twofl_get_pthermal_n
3550 
3551  subroutine twofl_get_pthermal_n_primitive(w,x,ixI^L,ixO^L,pth)
3553  integer, intent(in) :: ixI^L, ixO^L
3554  double precision, intent(in) :: w(ixI^S,1:nw)
3555  double precision, intent(in) :: x(ixI^S,1:ndim)
3556  double precision, intent(out) :: pth(ixI^S)
3557 
3558  if(phys_energy) then
3559  if(has_equi_pe_n0) then
3560  pth(ixo^s) = w(ixo^s,e_n_) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
3561  else
3562  pth(ixo^s) = w(ixo^s,e_n_)
3563  endif
3564  else
3565  call get_rhon_tot(w,x,ixi^l,ixo^l,pth)
3566  pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
3567  end if
3568  end subroutine twofl_get_pthermal_n_primitive
3569 
3570  !> Calculate v component
3571  subroutine twofl_get_v_n_idim(w,x,ixI^L,ixO^L,idim,v)
3573 
3574  integer, intent(in) :: ixi^l, ixo^l, idim
3575  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
3576  double precision, intent(out) :: v(ixi^s)
3577  double precision :: rhon(ixi^s)
3578 
3579  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
3580  v(ixo^s) = w(ixo^s, mom_n(idim)) / rhon(ixo^s)
3581 
3582  end subroutine twofl_get_v_n_idim
3583 
3584  subroutine internal_energy_add_source_n(qdt,ixI^L,ixO^L,wCT,w,x)
3586  use mod_geometry
3587 
3588  integer, intent(in) :: ixI^L, ixO^L
3589  double precision, intent(in) :: qdt
3590  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3591  double precision, intent(inout) :: w(ixI^S,1:nw)
3592  double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3593 
3594  call twofl_get_pthermal_n(wct,x,ixi^l,ixo^l,pth)
3595  call twofl_get_v_n(wct,x,ixi^l,ixi^l,v)
3596  call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-pth,w,x,e_n_)
3597 
3598  if(fix_small_values .and. .not. has_equi_pe_n0) then
3599  call twofl_handle_small_ei_n(w,x,ixi^l,ixo^l,e_n_,'internal_energy_add_source')
3600  end if
3601  end subroutine internal_energy_add_source_n
3602 
3603  !> Calculate v_c vector
3604  subroutine twofl_get_v_c(w,x,ixI^L,ixO^L,v)
3606 
3607  integer, intent(in) :: ixI^L, ixO^L
3608  double precision, intent(in) :: w(ixI^S,nw), x(ixI^S,1:ndim)
3609  double precision, intent(out) :: v(ixI^S,ndir)
3610  double precision :: rhoc(ixI^S)
3611  integer :: idir
3612 
3613  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3614  do idir=1,ndir
3615  v(ixo^s,idir) = w(ixo^s, mom_c(idir)) / rhoc(ixo^s)
3616  end do
3617 
3618  end subroutine twofl_get_v_c
3619 
3620  subroutine get_rhoc_tot(w,x,ixI^L,ixO^L,rhoc)
3622  integer, intent(in) :: ixi^l, ixo^l
3623  double precision, intent(in) :: w(ixi^s,1:nw), x(ixi^s,1:ndim)
3624  double precision, intent(out) :: rhoc(ixi^s)
3625  if(has_equi_rho_c0) then
3626  rhoc(ixo^s) = w(ixo^s,rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,b0i)
3627  else
3628  rhoc(ixo^s) = w(ixo^s,rho_c_)
3629  endif
3630 
3631  end subroutine get_rhoc_tot
3632 
3633  subroutine twofl_get_pthermal_c(w,x,ixI^L,ixO^L,pth)
3636  integer, intent(in) :: ixi^l, ixo^l
3637  double precision, intent(in) :: w(ixi^s,1:nw)
3638  double precision, intent(in) :: x(ixi^s,1:ndim)
3639  double precision, intent(out) :: pth(ixi^s)
3640  integer :: ix^d, iw
3641 
3642  if(phys_energy) then
3643  if(phys_internal_e) then
3644  pth(ixo^s)=gamma_1*w(ixo^s,e_c_)
3645  elseif(phys_total_energy) then
3646  pth(ixo^s)=gamma_1*(w(ixo^s,e_c_)&
3647  - twofl_kin_en_c(w,ixi^l,ixo^l)&
3648  - twofl_mag_en(w,ixi^l,ixo^l))
3649  else
3650  pth(ixo^s)=gamma_1*(w(ixo^s,e_c_)&
3651  - twofl_kin_en_c(w,ixi^l,ixo^l))
3652  end if
3653  if(has_equi_pe_c0) then
3654  pth(ixo^s) = pth(ixo^s) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
3655  endif
3656  else
3657  call get_rhoc_tot(w,x,ixi^l,ixo^l,pth)
3658  pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
3659  end if
3660 
3661  if (fix_small_values) then
3662  {do ix^db= ixo^lim^db\}
3663  if(pth(ix^d)<small_pressure) then
3664  pth(ix^d)=small_pressure
3665  end if
3666  {enddo^d&\}
3667  else if (check_small_values) then
3668  {do ix^db= ixo^lim^db\}
3669  if(pth(ix^d)<small_pressure) then
3670  write(*,*) "Error: small value of gas pressure",pth(ix^d),&
3671  " encountered when call twofl_get_pe_c1"
3672  write(*,*) "Iteration: ", it, " Time: ", global_time
3673  write(*,*) "Location: ", x(ix^d,:)
3674  write(*,*) "Cell number: ", ix^d
3675  do iw=1,nw
3676  write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
3677  end do
3678  ! use erroneous arithmetic operation to crash the run
3679  if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
3680  write(*,*) "Saving status at the previous time step"
3681  crash=.true.
3682  end if
3683  {enddo^d&\}
3684  end if
3685 
3686  end subroutine twofl_get_pthermal_c
3687 
3688  subroutine twofl_get_pthermal_c_primitive(w,x,ixI^L,ixO^L,pth)
3690  integer, intent(in) :: ixI^L, ixO^L
3691  double precision, intent(in) :: w(ixI^S,1:nw)
3692  double precision, intent(in) :: x(ixI^S,1:ndim)
3693  double precision, intent(out) :: pth(ixI^S)
3694 
3695  if(phys_energy) then
3696  if(has_equi_pe_c0) then
3697  pth(ixo^s) = w(ixo^s,e_c_) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
3698  else
3699  pth(ixo^s) = w(ixo^s,e_c_)
3700  endif
3701  else
3702  call get_rhoc_tot(w,x,ixi^l,ixo^l,pth)
3703  pth(ixo^s)=twofl_adiab*pth(ixo^s)**twofl_gamma
3704  end if
3705  end subroutine twofl_get_pthermal_c_primitive
3706 
3707  !> Calculate v_c component
3708  subroutine twofl_get_v_c_idim(w,x,ixI^L,ixO^L,idim,v)
3710 
3711  integer, intent(in) :: ixi^l, ixo^l, idim
3712  double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
3713  double precision, intent(out) :: v(ixi^s)
3714  double precision :: rhoc(ixi^s)
3715 
3716  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3717  v(ixo^s) = w(ixo^s, mom_c(idim)) / rhoc(ixo^s)
3718 
3719  end subroutine twofl_get_v_c_idim
3720 
3721  subroutine internal_energy_add_source_c(qdt,ixI^L,ixO^L,wCT,w,x,ie)
3723  use mod_geometry
3724 
3725  integer, intent(in) :: ixI^L, ixO^L,ie
3726  double precision, intent(in) :: qdt
3727  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3728  double precision, intent(inout) :: w(ixI^S,1:nw)
3729  double precision :: pth(ixI^S),v(ixI^S,1:ndir),divv(ixI^S)
3730 
3731  call twofl_get_pthermal_c(wct,x,ixi^l,ixo^l,pth)
3732  call twofl_get_v_c(wct,x,ixi^l,ixi^l,v)
3733  call add_geom_pdivv(qdt,ixi^l,ixo^l,v,-pth,w,x,ie)
3734  if(fix_small_values .and. .not. has_equi_pe_c0) then
3735  call twofl_handle_small_ei_c(w,x,ixi^l,ixo^l,ie,'internal_energy_add_source')
3736  end if
3737  end subroutine internal_energy_add_source_c
3738 
3739  !> handle small or negative internal energy
3740  subroutine twofl_handle_small_ei_c(w, x, ixI^L, ixO^L, ie, subname)
3742  use mod_small_values
3743  integer, intent(in) :: ixI^L,ixO^L, ie
3744  double precision, intent(inout) :: w(ixI^S,1:nw)
3745  double precision, intent(in) :: x(ixI^S,1:ndim)
3746  character(len=*), intent(in) :: subname
3747 
3748  integer :: idir
3749  logical :: flag(ixI^S,1:nw)
3750  double precision :: rhoc(ixI^S)
3751  double precision :: rhon(ixI^S)
3752 
3753  flag=.false.
3754  if(has_equi_pe_c0) then
3755  where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1<small_e)&
3756  flag(ixo^s,ie)=.true.
3757  else
3758  where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3759  endif
3760  if(any(flag(ixo^s,ie))) then
3761  select case (small_values_method)
3762  case ("replace")
3763  if(has_equi_pe_c0) then
3764  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3765  block%equi_vars(ixo^s,equi_pe_c0_,0)*inv_gamma_1
3766  else
3767  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3768  endif
3769  case ("average")
3770  call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
3771  case default
3772  ! small values error shows primitive variables
3773  ! to_primitive subroutine cannot be used as this error handling
3774  ! is also used in TC where e_to_ei is explicitly called
3775  w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
3776  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
3777  w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
3778  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3779  do idir = 1, ndir
3780  w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
3781  w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
3782  end do
3783  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
3784  end select
3785  end if
3786 
3787  end subroutine twofl_handle_small_ei_c
3788 
3789  !> handle small or negative internal energy
3790  subroutine twofl_handle_small_ei_n(w, x, ixI^L, ixO^L, ie, subname)
3792  use mod_small_values
3793  integer, intent(in) :: ixI^L,ixO^L, ie
3794  double precision, intent(inout) :: w(ixI^S,1:nw)
3795  double precision, intent(in) :: x(ixI^S,1:ndim)
3796  character(len=*), intent(in) :: subname
3797 
3798  integer :: idir
3799  logical :: flag(ixI^S,1:nw)
3800  double precision :: rhoc(ixI^S)
3801  double precision :: rhon(ixI^S)
3802 
3803  flag=.false.
3804  if(has_equi_pe_n0) then
3805  where(w(ixo^s,ie)+block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1<small_e)&
3806  flag(ixo^s,ie)=.true.
3807  else
3808  where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
3809  endif
3810  if(any(flag(ixo^s,ie))) then
3811  select case (small_values_method)
3812  case ("replace")
3813  if(has_equi_pe_n0) then
3814  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e - &
3815  block%equi_vars(ixo^s,equi_pe_n0_,0)*inv_gamma_1
3816  else
3817  where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
3818  endif
3819  case ("average")
3820  call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
3821  case default
3822  ! small values error shows primitive variables
3823  w(ixo^s,e_n_)=w(ixo^s,e_n_)*gamma_1
3824  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
3825  w(ixo^s,e_c_)=w(ixo^s,e_c_)*gamma_1
3826  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
3827  do idir = 1, ndir
3828  w(ixo^s, mom_n(idir)) = w(ixo^s, mom_n(idir))/rhon(ixo^s)
3829  w(ixo^s, mom_c(idir)) = w(ixo^s, mom_c(idir))/rhoc(ixo^s)
3830  end do
3831  call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
3832  end select
3833  end if
3834 
3835  end subroutine twofl_handle_small_ei_n
3836 
3837  !> Source terms after split off time-independent magnetic field
3838  subroutine add_source_b0split(qdt,ixI^L,ixO^L,wCT,w,x)
3840 
3841  integer, intent(in) :: ixI^L, ixO^L
3842  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3843  double precision, intent(inout) :: w(ixI^S,1:nw)
3844 
3845  double precision :: a(ixI^S,3), b(ixI^S,3), axb(ixI^S,3)
3846  integer :: idir
3847 
3848  a=0.d0
3849  b=0.d0
3850  ! for force-free field J0xB0 =0
3851  if(.not.b0field_forcefree) then
3852  ! store B0 magnetic field in b
3853  b(ixo^s,1:ndir)=block%B0(ixo^s,1:ndir,0)
3854 
3855  ! store J0 current in a
3856  do idir=7-2*ndir,3
3857  a(ixo^s,idir)=block%J0(ixo^s,idir)
3858  end do
3859  call cross_product(ixi^l,ixo^l,a,b,axb)
3860  axb(ixo^s,:)=axb(ixo^s,:)*qdt
3861  ! add J0xB0 source term in momentum equations
3862  w(ixo^s,mom_c(1:ndir))=w(ixo^s,mom_c(1:ndir))+axb(ixo^s,1:ndir)
3863  end if
3864 
3865  if(phys_total_energy) then
3866  a=0.d0
3867  ! for free-free field -(vxB0) dot J0 =0
3868  b(ixo^s,:)=wct(ixo^s,mag(:))
3869  ! store full magnetic field B0+B1 in b
3870  if(.not.b0field_forcefree) b(ixo^s,:)=b(ixo^s,:)+block%B0(ixo^s,:,0)
3871  ! store velocity in a
3872  do idir=1,ndir
3873  call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idir,a(ixi^s,idir))
3874  end do
3875  call cross_product(ixi^l,ixo^l,a,b,axb)
3876  axb(ixo^s,:)=axb(ixo^s,:)*qdt
3877  ! add -(vxB) dot J0 source term in energy equation
3878  do idir=7-2*ndir,3
3879  w(ixo^s,e_c_)=w(ixo^s,e_c_)-axb(ixo^s,idir)*block%J0(ixo^s,idir)
3880  end do
3881  end if
3882 
3883  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_B0')
3884 
3885  end subroutine add_source_b0split
3886 
3887  !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
3888  !> each direction, non-conservative. If the fourthorder precompiler flag is
3889  !> set, uses fourth order central difference for the laplacian. Then the
3890  !> stencil is 5 (2 neighbours).
3891  subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
3893  use mod_usr_methods
3894  use mod_geometry
3895 
3896  integer, intent(in) :: ixI^L, ixO^L
3897  double precision, intent(in) :: qdt
3898  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
3899  double precision, intent(inout) :: w(ixI^S,1:nw)
3900  integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
3901  integer :: lxO^L, kxO^L
3902 
3903  double precision :: tmp(ixI^S),tmp2(ixI^S)
3904 
3905  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
3906  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
3907  double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
3908 
3909  ! Calculating resistive sources involve one extra layer
3910  if (twofl_4th_order) then
3911  ixa^l=ixo^l^ladd2;
3912  else
3913  ixa^l=ixo^l^ladd1;
3914  end if
3915 
3916  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
3917  call mpistop("Error in add_source_res1: Non-conforming input limits")
3918 
3919  ! Calculate current density and idirmin
3920  call get_current(wct,ixi^l,ixo^l,idirmin,current)
3921 
3922  if (twofl_eta>zero)then
3923  eta(ixa^s)=twofl_eta
3924  gradeta(ixo^s,1:ndim)=zero
3925  else
3926  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
3927  ! assumes that eta is not function of current?
3928  do idim=1,ndim
3929  call gradient(eta,ixi^l,ixo^l,idim,tmp)
3930  gradeta(ixo^s,idim)=tmp(ixo^s)
3931  end do
3932  end if
3933 
3934  if(b0field) then
3935  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))+block%B0(ixi^s,1:ndir,0)
3936  else
3937  bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
3938  end if
3939 
3940  do idir=1,ndir
3941  ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
3942  if (twofl_4th_order) then
3943  tmp(ixo^s)=zero
3944  tmp2(ixi^s)=bf(ixi^s,idir)
3945  do idim=1,ndim
3946  lxo^l=ixo^l+2*kr(idim,^d);
3947  jxo^l=ixo^l+kr(idim,^d);
3948  hxo^l=ixo^l-kr(idim,^d);
3949  kxo^l=ixo^l-2*kr(idim,^d);
3950  tmp(ixo^s)=tmp(ixo^s)+&
3951  (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
3952  /(12.0d0 * dxlevel(idim)**2)
3953  end do
3954  else
3955  tmp(ixo^s)=zero
3956  tmp2(ixi^s)=bf(ixi^s,idir)
3957  do idim=1,ndim
3958  jxo^l=ixo^l+kr(idim,^d);
3959  hxo^l=ixo^l-kr(idim,^d);
3960  tmp(ixo^s)=tmp(ixo^s)+&
3961  (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
3962  end do
3963  end if
3964 
3965  ! Multiply by eta
3966  tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
3967 
3968  ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
3969  if (twofl_eta<zero)then
3970  do jdir=1,ndim; do kdir=idirmin,3
3971  if (lvc(idir,jdir,kdir)/=0)then
3972  if (lvc(idir,jdir,kdir)==1)then
3973  tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3974  else
3975  tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
3976  end if
3977  end if
3978  end do; end do
3979  end if
3980 
3981  ! Add sources related to eta*laplB-grad(eta) x J to B and e
3982  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
3983  if (phys_total_energy) then
3984  w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*tmp(ixo^s)*bf(ixo^s,idir)
3985  end if
3986  end do ! idir
3987 
3988  if (phys_energy) then
3989  ! de/dt+=eta*J**2
3990  w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
3991  end if
3992 
3993  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res1')
3994 
3995  end subroutine add_source_res1
3996 
3997  !> Add resistive source to w within ixO
3998  !> Uses 5 point stencil (2 neighbours) in each direction, conservative
3999  subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
4001  use mod_usr_methods
4002  use mod_geometry
4003 
4004  integer, intent(in) :: ixI^L, ixO^L
4005  double precision, intent(in) :: qdt
4006  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4007  double precision, intent(inout) :: w(ixI^S,1:nw)
4008 
4009  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4010  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
4011  double precision :: tmpvec(ixI^S,1:3),tmp(ixO^S)
4012  integer :: ixA^L,idir,idirmin,idirmin1
4013 
4014  ixa^l=ixo^l^ladd2;
4015 
4016  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4017  call mpistop("Error in add_source_res2: Non-conforming input limits")
4018 
4019  ixa^l=ixo^l^ladd1;
4020  ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
4021  ! Determine exact value of idirmin while doing the loop.
4022  call get_current(wct,ixi^l,ixa^l,idirmin,current)
4023 
4024  if (twofl_eta>zero)then
4025  eta(ixa^s)=twofl_eta
4026  else
4027  call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
4028  end if
4029 
4030  ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
4031  tmpvec(ixa^s,1:ndir)=zero
4032  do idir=idirmin,3
4033  tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
4034  end do
4035  curlj=0.d0
4036  call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
4037  if(stagger_grid.and.ndim==2.and.ndir==3) then
4038  ! if 2.5D
4039  w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
4040  else
4041  w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
4042  end if
4043 
4044  if(phys_energy) then
4045  if(phys_total_energy) then
4046  ! de/dt= +div(B x Jeta) = eta J^2 - B dot curl(eta J)
4047  ! de1/dt= eta J^2 - B1 dot curl(eta J)
4048  w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*(eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)-&
4049  sum(wct(ixo^s,mag(1:ndir))*curlj(ixo^s,1:ndir),dim=ndim+1))
4050  else
4051  ! add eta*J**2 source term in the internal energy equation
4052  w(ixo^s,e_c_)=w(ixo^s,e_c_)+qdt*eta(ixo^s)*sum(current(ixo^s,:)**2,dim=ndim+1)
4053  end if
4054 
4055  end if
4056 
4057  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_res2')
4058  end subroutine add_source_res2
4059 
4060  !> Add Hyper-resistive source to w within ixO
4061  !> Uses 9 point stencil (4 neighbours) in each direction.
4062  subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
4064  use mod_geometry
4065 
4066  integer, intent(in) :: ixI^L, ixO^L
4067  double precision, intent(in) :: qdt
4068  double precision, intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4069  double precision, intent(inout) :: w(ixI^S,1:nw)
4070  !.. local ..
4071  double precision :: current(ixI^S,7-2*ndir:3)
4072  double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
4073  integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
4074 
4075  ixa^l=ixo^l^ladd3;
4076  if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
4077  call mpistop("Error in add_source_hyperres: Non-conforming input limits")
4078 
4079  call get_current(wct,ixi^l,ixa^l,idirmin,current)
4080  tmpvec(ixa^s,1:ndir)=zero
4081  do jdir=idirmin,3
4082  tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
4083  end do
4084 
4085  ixa^l=ixo^l^ladd2;
4086  call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4087 
4088  ixa^l=ixo^l^ladd1;
4089  tmpvec(ixa^s,1:ndir)=zero
4090  call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
4091  ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*twofl_eta_hyper
4092 
4093  ixa^l=ixo^l;
4094  tmpvec2(ixa^s,1:ndir)=zero
4095  call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
4096 
4097  do idir=1,ndir
4098  w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
4099  end do
4100 
4101  if (phys_energy) then
4102  ! de/dt= +div(B x Ehyper)
4103  ixa^l=ixo^l^ladd1;
4104  tmpvec2(ixa^s,1:ndir)=zero
4105  do idir=1,ndir; do jdir=1,ndir; do kdir=idirmin,3
4106  tmpvec2(ixa^s,idir) = tmpvec(ixa^s,idir)&
4107  + lvc(idir,jdir,kdir)*wct(ixa^s,mag(jdir))*ehyper(ixa^s,kdir)
4108  end do; end do; end do
4109  tmp(ixo^s)=zero
4110  call divvector(tmpvec2,ixi^l,ixo^l,tmp)
4111  w(ixo^s,e_c_)=w(ixo^s,e_c_)+tmp(ixo^s)*qdt
4112  end if
4113 
4114  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_hyperres')
4115 
4116  end subroutine add_source_hyperres
4117 
4118  subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
4119  ! Add divB related sources to w within ixO
4120  ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
4121  ! giving the EGLM-MHD scheme
4123  use mod_geometry
4124 
4125  integer, intent(in) :: ixI^L, ixO^L
4126  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4127  double precision, intent(inout) :: w(ixI^S,1:nw)
4128  double precision:: divb(ixI^S)
4129  integer :: idim,idir
4130  double precision :: gradPsi(ixI^S)
4131 
4132  ! We calculate now div B
4133  call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
4134 
4135  ! dPsi/dt = - Ch^2/Cp^2 Psi
4136  if (twofl_glm_alpha < zero) then
4137  w(ixo^s,psi_) = abs(twofl_glm_alpha)*wct(ixo^s,psi_)
4138  else
4139  ! implicit update of Psi variable
4140  ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
4141  if(slab_uniform) then
4142  w(ixo^s,psi_) = dexp(-qdt*cmax_global*twofl_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
4143  else
4144  w(ixo^s,psi_) = dexp(-qdt*cmax_global*twofl_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
4145  end if
4146  end if
4147 
4148  ! gradient of Psi
4149  do idim=1,ndim
4150  select case(typegrad)
4151  case("central")
4152  call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
4153  case("limited")
4154  call gradients(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
4155  end select
4156  if (phys_total_energy) then
4157  ! e = e -qdt (b . grad(Psi))
4158  w(ixo^s,e_c_) = w(ixo^s,e_c_)-qdt*wct(ixo^s,mag(idim))*gradpsi(ixo^s)
4159  end if
4160  end do
4161 
4162  ! m = m - qdt b div b
4163  do idir=1,ndir
4164  w(ixo^s,mom_c(idir))=w(ixo^s,mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
4165  end do
4166 
4167  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_glm')
4168 
4169  end subroutine add_source_glm
4170 
4171  !> Add divB related sources to w within ixO corresponding to Powel
4172  subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
4174 
4175  integer, intent(in) :: ixI^L, ixO^L
4176  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4177  double precision, intent(inout) :: w(ixI^S,1:nw)
4178  double precision :: divb(ixI^S),v(ixI^S,1:ndir)
4179  integer :: idir
4180 
4181  ! We calculate now div B
4182  call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
4183 
4184  ! calculate velocity
4185  call twofl_get_v_c(wct,x,ixi^l,ixo^l,v)
4186 
4187  if (phys_total_energy) then
4188  ! e = e - qdt (v . b) * div b
4189  w(ixo^s,e_c_)=w(ixo^s,e_c_)-&
4190  qdt*sum(v(ixo^s,:)*wct(ixo^s,mag(:)),dim=ndim+1)*divb(ixo^s)
4191  end if
4192 
4193  ! b = b - qdt v * div b
4194  do idir=1,ndir
4195  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*v(ixo^s,idir)*divb(ixo^s)
4196  end do
4197 
4198  ! m = m - qdt b div b
4199  do idir=1,ndir
4200  w(ixo^s,mom_c(idir))=w(ixo^s,mom_c(idir))-qdt*twofl_mag_i_all(w,ixi^l,ixo^l,idir)*divb(ixo^s)
4201  end do
4202 
4203  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_powel')
4204 
4205  end subroutine add_source_powel
4206 
4207  subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
4208  ! Add divB related sources to w within ixO
4209  ! corresponding to Janhunen, just the term in the induction equation.
4211 
4212  integer, intent(in) :: ixI^L, ixO^L
4213  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4214  double precision, intent(inout) :: w(ixI^S,1:nw)
4215  double precision :: divb(ixI^S),vel(ixI^S)
4216  integer :: idir
4217 
4218  ! We calculate now div B
4219  call get_divb(wct,ixi^l,ixo^l,divb, twofl_divb_4thorder)
4220 
4221  ! b = b - qdt v * div b
4222  do idir=1,ndir
4223  call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idir,vel)
4224  w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*vel(ixo^s)*divb(ixo^s)
4225  end do
4226 
4227  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_janhunen')
4228 
4229  end subroutine add_source_janhunen
4230 
4231  subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
4232  ! Add Linde's divB related sources to wnew within ixO
4234  use mod_geometry
4235 
4236  integer, intent(in) :: ixI^L, ixO^L
4237  double precision, intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
4238  double precision, intent(inout) :: w(ixI^S,1:nw)
4239  integer :: idim, idir, ixp^L, i^D, iside
4240  double precision :: divb(ixI^S),graddivb(ixI^S)
4241  logical, dimension(-1:1^D&) :: leveljump
4242 
4243  ! Calculate div B
4244  ixp^l=ixo^l^ladd1;
4245  call get_divb(wct,ixi^l,ixp^l,divb, twofl_divb_4thorder)
4246 
4247  ! for AMR stability, retreat one cell layer from the boarders of level jump
4248  {do i^db=-1,1\}
4249  if(i^d==0|.and.) cycle
4250  if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
4251  leveljump(i^d)=.true.
4252  else
4253  leveljump(i^d)=.false.
4254  end if
4255  {end do\}
4256 
4257  ixp^l=ixo^l;
4258  do idim=1,ndim
4259  select case(idim)
4260  {case(^d)
4261  do iside=1,2
4262  i^dd=kr(^dd,^d)*(2*iside-3);
4263  if (leveljump(i^dd)) then
4264  if (iside==1) then
4265  ixpmin^d=ixomin^d-i^d
4266  else
4267  ixpmax^d=ixomax^d-i^d
4268  end if
4269  end if
4270  end do
4271  \}
4272  end select
4273  end do
4274 
4275  ! Add Linde's diffusive terms
4276  do idim=1,ndim
4277  ! Calculate grad_idim(divb)
4278  select case(typegrad)
4279  case("central")
4280  call gradient(divb,ixi^l,ixp^l,idim,graddivb)
4281  case("limited")
4282  call gradients(divb,ixi^l,ixp^l,idim,graddivb)
4283  end select
4284 
4285  ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
4286  if (slab_uniform) then
4287  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
4288  else
4289  graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
4290  /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
4291  end if
4292 
4293  w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
4294 
4295  if (typedivbdiff=='all' .and. phys_total_energy) then
4296  ! e += B_idim*eta*grad_idim(divb)
4297  w(ixp^s,e_c_)=w(ixp^s,e_c_)+wct(ixp^s,mag(idim))*graddivb(ixp^s)
4298  end if
4299  end do
4300 
4301  if (fix_small_values) call twofl_handle_small_values(.false.,w,x,ixi^l,ixo^l,'add_source_linde')
4302 
4303  end subroutine add_source_linde
4304 
4305 
4306  !> get dimensionless div B = |divB| * volume / area / |B|
4307  subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
4308 
4310 
4311  integer, intent(in) :: ixi^l, ixo^l
4312  double precision, intent(in) :: w(ixi^s,1:nw)
4313  double precision :: divb(ixi^s), dsurface(ixi^s)
4314 
4315  double precision :: invb(ixo^s)
4316  integer :: ixa^l,idims
4317 
4318  call get_divb(w,ixi^l,ixo^l,divb)
4319  invb(ixo^s)=sqrt(twofl_mag_en_all(w,ixi^l,ixo^l))
4320  where(invb(ixo^s)/=0.d0)
4321  invb(ixo^s)=1.d0/invb(ixo^s)
4322  end where
4323  if(slab_uniform) then
4324  divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
4325  else
4326  ixamin^d=ixomin^d-1;
4327  ixamax^d=ixomax^d-1;
4328  dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
4329  do idims=1,ndim
4330  ixa^l=ixo^l-kr(idims,^d);
4331  dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
4332  end do
4333  divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
4334  block%dvolume(ixo^s)/dsurface(ixo^s)
4335  end if
4336 
4337  end subroutine get_normalized_divb
4338 
4339  !> Calculate idirmin and the idirmin:3 components of the common current array
4340  !> make sure that dxlevel(^D) is set correctly.
4341  subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
4343  use mod_geometry
4344 
4345  integer, intent(in) :: ixo^l, ixi^l
4346  double precision, intent(in) :: w(ixi^s,1:nw)
4347  integer, intent(out) :: idirmin
4348  integer :: idir, idirmin0
4349 
4350  ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
4351  double precision :: current(ixi^s,7-2*ndir:3),bvec(ixi^s,1:ndir)
4352 
4353  idirmin0 = 7-2*ndir
4354 
4355  bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
4356 
4357  call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
4358 
4359  if(b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
4360  block%J0(ixo^s,idirmin0:3)
4361 
4362  end subroutine get_current
4363 
4364  ! copied from gravity
4365  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
4366  subroutine gravity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
4367  energy,qsourcesplit,active)
4369  use mod_usr_methods
4370 
4371  integer, intent(in) :: ixI^L, ixO^L
4372  double precision, intent(in) :: qdt, x(ixI^S,1:ndim)
4373  double precision, intent(in) :: wCT(ixI^S,1:nw)
4374  double precision, intent(inout) :: w(ixI^S,1:nw)
4375  logical, intent(in) :: energy,qsourcesplit
4376  logical, intent(inout) :: active
4377  double precision :: vel(ixI^S)
4378  integer :: idim
4379 
4380  double precision :: gravity_field(ixI^S,ndim)
4381 
4382  if(qsourcesplit .eqv. grav_split) then
4383  active = .true.
4384 
4385  if (.not. associated(usr_gravity)) then
4386  write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
4387  write(*,*) "like the phys_gravity in mod_usr_methods.t"
4388  call mpistop("gravity_add_source: usr_gravity not defined")
4389  else
4390  call usr_gravity(ixi^l,ixo^l,wct,x,gravity_field)
4391  end if
4392 
4393  do idim = 1, ndim
4394  w(ixo^s,mom_n(idim)) = w(ixo^s,mom_n(idim)) &
4395  + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,rho_n_)
4396  w(ixo^s,mom_c(idim)) = w(ixo^s,mom_c(idim)) &
4397  + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,rho_c_)
4398  if(energy) then
4399 #if !defined(E_RM_W0) || E_RM_W0 == 1
4400  call twofl_get_v_n_idim(wct,x,ixi^l,ixo^l,idim,vel)
4401  w(ixo^s,e_n_)=w(ixo^s,e_n_) &
4402  + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,rho_n_)
4403  call twofl_get_v_c_idim(wct,x,ixi^l,ixo^l,idim,vel)
4404  w(ixo^s,e_c_)=w(ixo^s,e_c_) &
4405  + qdt * gravity_field(ixo^s,idim) * vel(ixo^s) * wct(ixo^s,rho_c_)
4406 #else
4407  w(ixo^s,e_n_)=w(ixo^s,e_n_) &
4408  + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,mom_n(idim))
4409  w(ixo^s,e_c_)=w(ixo^s,e_c_) &
4410  + qdt * gravity_field(ixo^s,idim) * wct(ixo^s,mom_c(idim))
4411 #endif
4412 
4413 
4414  end if
4415  end do
4416  end if
4417 
4418  end subroutine gravity_add_source
4419 
4420  subroutine gravity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4422  use mod_usr_methods
4423 
4424  integer, intent(in) :: ixI^L, ixO^L
4425  double precision, intent(in) :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
4426  double precision, intent(inout) :: dtnew
4427 
4428  double precision :: dxinv(1:ndim), max_grav
4429  integer :: idim
4430 
4431  double precision :: gravity_field(ixI^S,ndim)
4432 
4433  ^d&dxinv(^d)=one/dx^d;
4434 
4435  if(.not. associated(usr_gravity)) then
4436  write(*,*) "mod_usr.t: please point usr_gravity to a subroutine"
4437  write(*,*) "like the phys_gravity in mod_usr_methods.t"
4438  call mpistop("gravity_get_dt: usr_gravity not defined")
4439  else
4440  call usr_gravity(ixi^l,ixo^l,w,x,gravity_field)
4441  end if
4442 
4443  do idim = 1, ndim
4444  max_grav = maxval(abs(gravity_field(ixo^s,idim)))
4445  max_grav = max(max_grav, epsilon(1.0d0))
4446  dtnew = min(dtnew, 1.0d0 / sqrt(max_grav * dxinv(idim)))
4447  end do
4448 
4449  end subroutine gravity_get_dt
4450 
4451  !> If resistivity is not zero, check diffusion time limit for dt
4452  subroutine twofl_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4454  use mod_usr_methods
4456  !use mod_viscosity, only: viscosity_get_dt
4457  !use mod_gravity, only: gravity_get_dt
4458 
4459  integer, intent(in) :: ixI^L, ixO^L
4460  double precision, intent(inout) :: dtnew
4461  double precision, intent(in) :: dx^D
4462  double precision, intent(in) :: w(ixI^S,1:nw)
4463  double precision, intent(in) :: x(ixI^S,1:ndim)
4464 
4465  integer :: idirmin,idim
4466  double precision :: dxarr(ndim)
4467  double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
4468 
4469  dtnew = bigdouble
4470 
4471  ^d&dxarr(^d)=dx^d;
4472  if (twofl_eta>zero)then
4473  dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/twofl_eta
4474  else if (twofl_eta<zero)then
4475  call get_current(w,ixi^l,ixo^l,idirmin,current)
4476  call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
4477  dtnew=bigdouble
4478  do idim=1,ndim
4479  if(slab_uniform) then
4480  dtnew=min(dtnew,&
4481  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
4482  else
4483  dtnew=min(dtnew,&
4484  dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
4485  end if
4486  end do
4487  end if
4488 
4489  if(twofl_eta_hyper>zero) then
4490  if(slab_uniform) then
4491  dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/twofl_eta_hyper,dtnew)
4492  else
4493  dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/twofl_eta_hyper,dtnew)
4494  end if
4495  end if
4496 
4497  ! the timestep related to coll terms: 1/(rho_n rho_c alpha)
4498  if(dtcollpar>0d0 .and. has_collisions()) then
4499  call coll_get_dt(w,x,ixi^l,ixo^l,dtnew)
4500  endif
4501 
4502  if(twofl_radiative_cooling_c) then
4503  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl_c)
4504  end if
4505  if(twofl_radiative_cooling_n) then
4506  call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl_n)
4507  end if
4508 !
4509 ! if(twofl_viscosity) then
4510 ! call viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
4511 ! end if
4512 !
4513  if(twofl_gravity) then
4514  call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
4515  end if
4516  if(twofl_hyperdiffusivity) then
4517  call hyperdiffusivity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
4518  end if
4519 
4520 
4521  end subroutine twofl_get_dt
4522 
4523  pure function has_collisions() result(res)
4524  logical :: res
4525  res = .not. twofl_alpha_coll_constant .or. twofl_alpha_coll >0d0
4526  end function has_collisions
4527 
4528  subroutine coll_get_dt(w,x,ixI^L,ixO^L,dtnew)
4530  integer, intent(in) :: ixi^l, ixo^l
4531  double precision, intent(in) :: w(ixi^s,1:nw)
4532  double precision, intent(in) :: x(ixi^s,1:ndim)
4533  double precision, intent(inout) :: dtnew
4534 
4535  double precision :: rhon(ixi^s), rhoc(ixi^s), alpha(ixi^s)
4536  double precision, allocatable :: gamma_rec(:^d&), gamma_ion(:^D&)
4537  double precision :: max_coll_rate
4538 
4539  call get_rhon_tot(w,x,ixi^l,ixo^l,rhon)
4540  call get_rhoc_tot(w,x,ixi^l,ixo^l,rhoc)
4541 
4542  call get_alpha_coll(ixi^l, ixo^l, w, x, alpha)
4543  max_coll_rate = maxval(alpha(ixo^s) * max(rhon(ixo^s), rhoc(ixo^s)))
4544 
4545  if(twofl_coll_inc_ionrec) then
4546  allocate(gamma_ion(ixi^s), gamma_rec(ixi^s))
4547  call get_gamma_ion_rec(ixi^l, ixo^l, w, x, gamma_rec, gamma_ion)
4548  max_coll_rate=max(max_coll_rate, maxval(gamma_ion(ixo^s)), maxval(gamma_rec(ixo^s)))
4549  deallocate(gamma_ion, gamma_rec)
4550  endif
4551  dtnew = min(dtcollpar/max_coll_rate, dtnew)
4552 
4553  end subroutine coll_get_dt
4554 
4555  ! Add geometrical source terms to w
4556  subroutine twofl_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
4558  use mod_geometry
4559 
4560  integer, intent(in) :: ixI^L, ixO^L
4561  double precision, intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
4562  double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
4563 
4564  integer :: iw,idir, h1x^L{^NOONED, h2x^L}
4565  double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),rho(ixI^S)
4566 
4567  integer :: mr_,mphi_ ! Polar var. names
4568  integer :: br_,bphi_
4569 
4570  ! charges
4571 
4572  mr_=mom_c(1); mphi_=mom_c(1)-1+phi_ ! Polar var. names
4573  br_=mag(1); bphi_=mag(1)-1+phi_
4574  call get_rhoc_tot(wct,x,ixi^l,ixo^l,rho)
4575 
4576  select case (coordinate)
4577  case (cylindrical)
4578  call twofl_get_p_c_total(wct,x,ixi^l,ixo^l,tmp)
4579 
4580  if(phi_>0) then
4581  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*(tmp(ixo^s)-&
4582  wct(ixo^s,bphi_)**2+wct(ixo^s,mphi_)**2/rho(ixo^s))
4583  w(ixo^s,mphi_)=w(ixo^s,mphi_)+qdt/x(ixo^s,1)*(&
4584  -wct(ixo^s,mphi_)*wct(ixo^s,mr_)/rho(ixo^s) &
4585  +wct(ixo^s,bphi_)*wct(ixo^s,br_))
4586  if(.not.stagger_grid) then
4587  w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt/x(ixo^s,1)*&
4588  (wct(ixo^s,bphi_)*wct(ixo^s,mr_) &
4589  -wct(ixo^s,br_)*wct(ixo^s,mphi_)) &
4590  /rho(ixo^s)
4591  end if
4592  else
4593  w(ixo^s,mr_)=w(ixo^s,mr_)+qdt/x(ixo^s,1)*tmp(ixo^s)
4594  end if
4595  if(twofl_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
4596  case (spherical)
4597  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
4598  call twofl_get_p_c_total(wct,x,ixi^l,ixo^l,tmp1)
4599  tmp(ixo^s)=tmp1(ixo^s)
4600  if(b0field) then
4601  tmp2(ixo^s)=sum(block%B0(ixo^s,:,0)*wct(ixo^s,mag(:)),dim=ndim+1)
4602  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4603  end if
4604  ! m1
4605  tmp(ixo^s)=tmp(ixo^s)*x(ixo^s,1) &
4606  *(block%surfaceC(ixo^s,1)-block%surfaceC(h1x^s,1))/block%dvolume(ixo^s)
4607  if(ndir>1) then
4608  do idir=2,ndir
4609  tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom_c(idir))**2/rho(ixo^s)-wct(ixo^s,mag(idir))**2
4610  if(b0field) tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,idir,0)*wct(ixo^s,mag(idir))
4611  end do
4612  end if
4613  w(ixo^s,mom_c(1))=w(ixo^s,mom_c(1))+qdt*tmp(ixo^s)/x(ixo^s,1)
4614  ! b1
4615  if(twofl_glm) then
4616  w(ixo^s,mag(1))=w(ixo^s,mag(1))+qdt/x(ixo^s,1)*2.0d0*wct(ixo^s,psi_)
4617  end if
4618 
4619  {^nooned
4620  ! m2
4621  tmp(ixo^s)=tmp1(ixo^s)
4622  if(b0field) then
4623  tmp(ixo^s)=tmp(ixo^s)+tmp2(ixo^s)
4624  end if
4625  ! This will make hydrostatic p=const an exact solution
4626  w(ixo^s,mom_c(2))=w(ixo^s,mom_c(2))+qdt*tmp(ixo^s) &
4627  *(block%surfaceC(ixo^s,2)-block%surfaceC(h2x^s,2)) &
4628  /block%dvolume(ixo^s)
4629  tmp(ixo^s)=-(wct(ixo^s,mom_c(1))*wct(ixo^s,mom_c(2))/rho(ixo^s) &
4630  -wct(ixo^s,mag(1))*wct(ixo^s,mag(2)))
4631  if (b0field) then
4632  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(2)) &
4633  +wct(ixo^s,mag(1))*block%B0(ixo^s,2,0)
4634  end if
4635  if(ndir==3) then
4636  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(3))**2/rho(ixo^s) &
4637  -wct(ixo^s,mag(3))**2)*dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4638  if (b0field) then
4639  tmp(ixo^s)=tmp(ixo^s)-2.0d0*block%B0(ixo^s,3,0)*wct(ixo^s,mag(3))&
4640  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2))
4641  end if
4642  end if
4643  w(ixo^s,mom_c(2))=w(ixo^s,mom_c(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4644  ! b2
4645  if(.not.stagger_grid) then
4646  tmp(ixo^s)=(wct(ixo^s,mom_c(1))*wct(ixo^s,mag(2)) &
4647  -wct(ixo^s,mom_c(2))*wct(ixo^s,mag(1)))/rho(ixo^s)
4648  if(b0field) then
4649  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(1))*block%B0(ixo^s,2,0) &
4650  -wct(ixo^s,mom_c(2))*block%B0(ixo^s,1,0))/rho(ixo^s)
4651  end if
4652  if(twofl_glm) then
4653  tmp(ixo^s)=tmp(ixo^s) &
4654  + dcos(x(ixo^s,2))/dsin(x(ixo^s,2))*wct(ixo^s,psi_)
4655  end if
4656  w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
4657  end if
4658  }
4659 
4660  if(ndir==3) then
4661  ! m3
4662  tmp(ixo^s)=-(wct(ixo^s,mom_c(3))*wct(ixo^s,mom_c(1))/rho(ixo^s) &
4663  -wct(ixo^s,mag(3))*wct(ixo^s,mag(1))) {^nooned &
4664  -(wct(ixo^s,mom_c(2))*wct(ixo^s,mom_c(3))/rho(ixo^s) &
4665  -wct(ixo^s,mag(2))*wct(ixo^s,mag(3))) &
4666  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4667  if (b0field) then
4668  tmp(ixo^s)=tmp(ixo^s)+block%B0(ixo^s,1,0)*wct(ixo^s,mag(3)) &
4669  +wct(ixo^s,mag(1))*block%B0(ixo^s,3,0) {^nooned &
4670  +(block%B0(ixo^s,2,0)*wct(ixo^s,mag(3)) &
4671  +wct(ixo^s,mag(2))*block%B0(ixo^s,3,0)) &
4672  *dcos(x(ixo^s,2))/dsin(x(ixo^s,2)) }
4673  end if
4674  w(ixo^s,mom_c(3))=w(ixo^s,mom_c(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4675  ! b3
4676  if(.not.stagger_grid) then
4677  tmp(ixo^s)=(wct(ixo^s,mom_c(1))*wct(ixo^s,mag(3)) &
4678  -wct(ixo^s,mom_c(3))*wct(ixo^s,mag(1)))/rho(ixo^s) {^nooned &
4679  -(wct(ixo^s,mom_c(3))*wct(ixo^s,mag(2)) &
4680  -wct(ixo^s,mom_c(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
4681  /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4682  if (b0field) then
4683  tmp(ixo^s)=tmp(ixo^s)+(wct(ixo^s,mom_c(1))*block%B0(ixo^s,3,0) &
4684  -wct(ixo^s,mom_c(3))*block%B0(ixo^s,1,0))/rho(ixo^s){^nooned &
4685  -(wct(ixo^s,mom_c(3))*block%B0(ixo^s,2,0) &
4686  -wct(ixo^s,mom_c(2))*block%B0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
4687  /(rho(ixo^s)*dsin(x(ixo^s,2))) }
4688  end if
4689  w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
4690  end if
4691  end if
4692  end select
4693 
4694  ! neutrals
4695  !TODO no dust: see and implement them from hd/mod_hd_phys !
4696  !uncomment cartesian expansion
4697  call get_rhon_tot(wct,x,ixi^l,ixo^l,rho)
4698  call twofl_get_pthermal_n(wct, x, ixi^l, ixo^l, tmp1)
4699 
4700  select case (coordinate)
4701 ! case(Cartesian_expansion)
4702 ! !the user provides the functions of exp_factor and del_exp_factor
4703 ! if(associated(usr_set_surface)) call usr_set_surface(ixI^L,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
4704 ! tmp(ixO^S) = tmp1(ixO^S)*del_exp_factor(ixO^S)/exp_factor(ixO^S)
4705 ! w(ixO^S,mom(1)) = w(ixO^S,mom(1)) + qdt*tmp(ixO^S)
4706 
4707  case (cylindrical)
4708  mr_ = mom_n(r_)
4709  if (phi_ > 0) then
4710  where (rho(ixo^s) > 0d0)
4711  tmp(ixo^s) = tmp1(ixo^s) + wct(ixo^s, mphi_)**2 / rho(ixo^s)
4712  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, r_)
4713  end where
4714  ! s[mphi]=(-mphi*mr/rho)/radius
4715  where (rho(ixo^s) > 0d0)
4716  tmp(ixo^s) = -wct(ixo^s, mphi_) * wct(ixo^s, mr_) / rho(ixo^s)
4717  w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt * tmp(ixo^s) / x(ixo^s, r_)
4718  end where
4719  else
4720  ! s[mr]=2pthermal/radius
4721  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp1(ixo^s) / x(ixo^s, r_)
4722  end if
4723  case (spherical)
4724  if(phi_>0) mphi_ = mom_n(phi_)
4725  h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
4726  ! s[mr]=((mtheta**2+mphi**2)/rho+2*p)/r
4727  tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4728  *(block%surfaceC(ixo^s, 1) - block%surfaceC(h1x^s, 1)) &
4729  /block%dvolume(ixo^s)
4730  if (ndir > 1) then
4731  do idir = 2, ndir
4732  tmp(ixo^s) = tmp(ixo^s) + wct(ixo^s, mom_n(idir))**2 / rho(ixo^s)
4733  end do
4734  end if
4735  w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4736 
4737  {^nooned
4738  ! s[mtheta]=-(mr*mtheta/rho)/r+cot(theta)*(mphi**2/rho+p)/r
4739  tmp(ixo^s) = tmp1(ixo^s) * x(ixo^s, 1) &
4740  * (block%surfaceC(ixo^s, 2) - block%surfaceC(h2x^s, 2)) &
4741  / block%dvolume(ixo^s)
4742  if (ndir == 3) then
4743  tmp(ixo^s) = tmp(ixo^s) + (wct(ixo^s, mom_n(3))**2 / rho(ixo^s)) / tan(x(ixo^s, 2))
4744  end if
4745  tmp(ixo^s) = tmp(ixo^s) - (wct(ixo^s, mom_n(2)) * wct(ixo^s, mr_)) / rho(ixo^s)
4746  w(ixo^s, mom_n(2)) = w(ixo^s, mom_n(2)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4747 
4748  if (ndir == 3) then
4749  ! s[mphi]=-(mphi*mr/rho)/r-cot(theta)*(mtheta*mphi/rho)/r
4750  tmp(ixo^s) = -(wct(ixo^s, mom_n(3)) * wct(ixo^s, mr_)) / rho(ixo^s)&
4751  - (wct(ixo^s, mom_n(2)) * wct(ixo^s, mom_n(3))) / rho(ixo^s) / tan(x(ixo^s, 2))
4752  w(ixo^s, mom_n(3)) = w(ixo^s, mom_n(3)) + qdt * tmp(ixo^s) / x(ixo^s, 1)
4753  end if
4754  }
4755  end select
4756 
4757  contains
4758  subroutine twofl_get_p_c_total(w,x,ixI^L,ixO^L,p)
4760 
4761  integer, intent(in) :: ixI^L, ixO^L
4762  double precision, intent(in) :: w(ixI^S,nw)
4763  double precision, intent(in) :: x(ixI^S,1:ndim)
4764  double precision, intent(out) :: p(ixI^S)
4765 
4766  call twofl_get_pthermal_c(w,x,ixi^l,ixo^l,p)
4767 
4768  p(ixo^s) = p(ixo^s) + 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4769 
4770  end subroutine twofl_get_p_c_total
4771 
4772  end subroutine twofl_add_source_geom
4773 
4774  subroutine twofl_get_temp_c_pert_from_etot(w, x, ixI^L, ixO^L, res)
4776  integer, intent(in) :: ixI^L, ixO^L
4777  double precision, intent(in) :: w(ixI^S, 1:nw)
4778  double precision, intent(in) :: x(ixI^S, 1:ndim)
4779  double precision, intent(out):: res(ixI^S)
4780 
4781  ! store pe1 in res
4782  res(ixo^s)=(gamma_1*(w(ixo^s,e_c_)&
4783  - twofl_kin_en_c(w,ixi^l,ixo^l)&
4784  - twofl_mag_en(w,ixi^l,ixo^l)))
4785  if(has_equi_pe_c0) then
4786  res(ixo^s) = res(ixo^s) + block%equi_vars(ixo^s,equi_pe_c0_,b0i)
4787  if(has_equi_rho_c0) then
4788  res(ixo^s) = res(ixo^s)/(rc * (w(ixo^s,rho_c_)+ block%equi_vars(ixo^s,equi_rho_c0_,b0i))) - &
4789  block%equi_vars(ixo^s,equi_pe_c0_,b0i)/(rc * block%equi_vars(ixo^s,equi_rho_c0_,b0i))
4790  else
4791  ! infinite equi temperature with p0 and 0 density
4792  res(ixo^s) = 0d0
4793  endif
4794  else
4795  res(ixo^s) = res(ixo^s)/(rc * w(ixo^s,rho_c_))
4796  endif
4797 
4798  end subroutine twofl_get_temp_c_pert_from_etot
4799 
4800  !> Compute 2 times total magnetic energy
4801  function twofl_mag_en_all(w, ixI^L, ixO^L) result(mge)
4803  integer, intent(in) :: ixi^l, ixo^l
4804  double precision, intent(in) :: w(ixi^s, nw)
4805  double precision :: mge(ixo^s)
4806 
4807  if (b0field) then
4808  mge(ixo^s) = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
4809  else
4810  mge(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4811  end if
4812  end function twofl_mag_en_all
4813 
4814  !> Compute full magnetic field by direction
4815  function twofl_mag_i_all(w, ixI^L, ixO^L,idir) result(mgf)
4817  integer, intent(in) :: ixi^l, ixo^l, idir
4818  double precision, intent(in) :: w(ixi^s, nw)
4819  double precision :: mgf(ixo^s)
4820 
4821  if (b0field) then
4822  mgf(ixo^s) = w(ixo^s, mag(idir))+block%B0(ixo^s,idir,b0i)
4823  else
4824  mgf(ixo^s) = w(ixo^s, mag(idir))
4825  end if
4826  end function twofl_mag_i_all
4827 
4828  !> Compute evolving magnetic energy
4829  function twofl_mag_en(w, ixI^L, ixO^L) result(mge)
4830  use mod_global_parameters, only: nw, ndim
4831  integer, intent(in) :: ixi^l, ixo^l
4832  double precision, intent(in) :: w(ixi^s, nw)
4833  double precision :: mge(ixo^s)
4834 
4835  mge(ixo^s) = 0.5d0 * sum(w(ixo^s, mag(:))**2, dim=ndim+1)
4836  end function twofl_mag_en
4837 
4838  !> compute kinetic energy of neutrals
4839  function twofl_kin_en_n(w, ixI^L, ixO^L) result(ke)
4840  use mod_global_parameters, only: nw, ndim,block
4841  integer, intent(in) :: ixi^l, ixo^l
4842  double precision, intent(in) :: w(ixi^s, nw)
4843  double precision :: ke(ixo^s)
4844 
4845  if(has_equi_rho_n0) then
4846  ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_n(:))**2, dim=ndim+1) / (w(ixo^s, rho_n_) + block%equi_vars(ixo^s,equi_rho_n0_,0))
4847  else
4848  ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_n(:))**2, dim=ndim+1) / w(ixo^s, rho_n_)
4849  endif
4850 
4851  end function twofl_kin_en_n
4852 
4853  subroutine twofl_get_temp_n_pert_from_etot(w, x, ixI^L, ixO^L, res)
4855  integer, intent(in) :: ixI^L, ixO^L
4856  double precision, intent(in) :: w(ixI^S, 1:nw)
4857  double precision, intent(in) :: x(ixI^S, 1:ndim)
4858  double precision, intent(out):: res(ixI^S)
4859 
4860  ! store pe1 in res
4861  res(ixo^s)=(gamma_1*(w(ixo^s,e_c_)- twofl_kin_en_c(w,ixi^l,ixo^l)))
4862  if(has_equi_pe_n0) then
4863  res(ixo^s) = res(ixo^s) + block%equi_vars(ixo^s,equi_pe_n0_,b0i)
4864  if(has_equi_rho_n0) then
4865  res(ixo^s) = res(ixo^s)/(rn * (w(ixo^s,rho_n_)+ block%equi_vars(ixo^s,equi_rho_n0_,b0i))) - &
4866  block%equi_vars(ixo^s,equi_pe_n0_,b0i)/(rn * block%equi_vars(ixo^s,equi_rho_n0_,b0i))
4867  else
4868  ! infinite equi temperature with p0 and 0 density
4869  res(ixo^s) = 0d0
4870  endif
4871  else
4872  res(ixo^s) = res(ixo^s)/(rn * w(ixo^s,rho_n_))
4873  endif
4874 
4875  end subroutine twofl_get_temp_n_pert_from_etot
4876 
4877  !> compute kinetic energy of charges
4878  !> w are conserved variables
4879  function twofl_kin_en_c(w, ixI^L, ixO^L) result(ke)
4880  use mod_global_parameters, only: nw, ndim,block
4881  integer, intent(in) :: ixi^l, ixo^l
4882  double precision, intent(in) :: w(ixi^s, nw)
4883  double precision :: ke(ixo^s)
4884 
4885  if(has_equi_rho_c0) then
4886  ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_c(:))**2, dim=ndim+1) / (w(ixo^s, rho_c_) + block%equi_vars(ixo^s,equi_rho_c0_,0))
4887  else
4888  ke(ixo^s) = 0.5d0 * sum(w(ixo^s, mom_c(:))**2, dim=ndim+1) / w(ixo^s, rho_c_)
4889  endif
4890  end function twofl_kin_en_c
4891 
4892  subroutine twofl_getv_hall(w,x,ixI^L,ixO^L,vHall)
4894 
4895  integer, intent(in) :: ixI^L, ixO^L
4896  double precision, intent(in) :: w(ixI^S,nw)
4897  double precision, intent(in) :: x(ixI^S,1:ndim)
4898  double precision, intent(inout) :: vHall(ixI^S,1:3)
4899 
4900  integer :: idir, idirmin
4901  double precision :: current(ixI^S,7-2*ndir:3)
4902  double precision :: rho(ixI^S)
4903 
4904  call get_rhoc_tot(w,x,ixi^l,ixo^l,rho)
4905  ! Calculate current density and idirmin
4906  call get_current(w,ixi^l,ixo^l,idirmin,current)
4907  vhall(ixo^s,1:3) = zero
4908  vhall(ixo^s,idirmin:3) = - twofl_etah*current(ixo^s,idirmin:3)
4909  do idir = idirmin, 3
4910  vhall(ixo^s,idir) = vhall(ixo^s,idir)/rho(ixo^s)
4911  end do
4912 
4913  end subroutine twofl_getv_hall
4914 
4915 ! the following not used
4916 ! subroutine twofl_getdt_Hall(w,x,ixI^L,ixO^L,dx^D,dthall)
4917 ! use mod_global_parameters
4918 !
4919 ! integer, intent(in) :: ixI^L, ixO^L
4920 ! double precision, intent(in) :: dx^D
4921 ! double precision, intent(in) :: w(ixI^S,1:nw)
4922 ! double precision, intent(in) :: x(ixI^S,1:ndim)
4923 ! double precision, intent(out) :: dthall
4924 ! !.. local ..
4925 ! double precision :: dxarr(ndim)
4926 ! double precision :: bmag(ixI^S)
4927 !
4928 ! dthall=bigdouble
4929 !
4930 ! ! because we have that in cmax now:
4931 ! return
4932 !
4933 ! ^D&dxarr(^D)=dx^D;
4934 !
4935 ! if (.not. B0field) then
4936 ! bmag(ixO^S)=sqrt(sum(w(ixO^S,mag(:))**2, dim=ndim+1))
4937 ! bmag(ixO^S)=sqrt(sum((w(ixO^S,mag(:)) + block%B0(ixO^S,1:ndir,b0i))**2))
4938 ! end if
4939 !
4940 ! if(slab_uniform) then
4941 ! dthall=dtdiffpar*minval(dxarr(1:ndim))**2.0d0/(twofl_etah*maxval(bmag(ixO^S)/w(ixO^S,rho_c_)))
4942 ! else
4943 ! dthall=dtdiffpar*minval(block%ds(ixO^S,1:ndim))**2.0d0/(twofl_etah*maxval(bmag(ixO^S)/w(ixO^S,rho_c_)))
4944 ! end if
4945 !
4946 ! end subroutine twofl_getdt_Hall
4947 
4948  subroutine twofl_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
4950  use mod_usr_methods
4951  integer, intent(in) :: ixI^L, ixO^L, idir
4952  double precision, intent(in) :: qt
4953  double precision, intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
4954  double precision, intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
4955  type(state) :: s
4956  double precision :: dB(ixI^S), dPsi(ixI^S)
4957 
4958  if(stagger_grid) then
4959  wlc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4960  wrc(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4961  wlp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4962  wrp(ixo^s,mag(idir))=s%ws(ixo^s,idir)
4963  else
4964  ! Solve the Riemann problem for the linear 2x2 system for normal
4965  ! B-field and GLM_Psi according to Dedner 2002:
4966  ! This implements eq. (42) in Dedner et al. 2002 JcP 175
4967  ! Gives the Riemann solution on the interface
4968  ! for the normal B component and Psi in the GLM-MHD system.
4969  ! 23/04/2013 Oliver Porth
4970  db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
4971  dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
4972 
4973  wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
4974  - 0.5d0/cmax_global * dpsi(ixo^s)
4975  wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
4976  - 0.5d0*cmax_global * db(ixo^s)
4977 
4978  wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4979  wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
4980 
4981  if(phys_total_energy) then
4982  wrc(ixo^s,e_c_)=wrc(ixo^s,e_c_)-half*wrc(ixo^s,mag(idir))**2
4983  wlc(ixo^s,e_c_)=wlc(ixo^s,e_c_)-half*wlc(ixo^s,mag(idir))**2
4984  end if
4985  wrc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4986  wrc(ixo^s,psi_) = wlp(ixo^s,psi_)
4987  wlc(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
4988  wlc(ixo^s,psi_) = wlp(ixo^s,psi_)
4989  ! modify total energy according to the change of magnetic field
4990  if(phys_total_energy) then
4991  wrc(ixo^s,e_c_)=wrc(ixo^s,e_c_)+half*wrc(ixo^s,mag(idir))**2
4992  wlc(ixo^s,e_c_)=wlc(ixo^s,e_c_)+half*wlc(ixo^s,mag(idir))**2
4993  end if
4994  end if
4995 
4996  if(associated(usr_set_wlr)) call usr_set_wlr(ixi^l,ixo^l,qt,wlc,wrc,wlp,wrp,s,idir)
4997 
4998  end subroutine twofl_modify_wlr
4999 
5000  subroutine twofl_boundary_adjust(igrid,psb)
5002  integer, intent(in) :: igrid
5003  type(state), target :: psb(max_blocks)
5004 
5005  integer :: iB, idims, iside, ixO^L, i^D
5006 
5007  block=>ps(igrid)
5008  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
5009  do idims=1,ndim
5010  ! to avoid using as yet unknown corner info in more than 1D, we
5011  ! fill only interior mesh ranges of the ghost cell ranges at first,
5012  ! and progressively enlarge the ranges to include corners later
5013  do iside=1,2
5014  i^d=kr(^d,idims)*(2*iside-3);
5015  if (neighbor_type(i^d,igrid)/=1) cycle
5016  ib=(idims-1)*2+iside
5017  if(.not.boundary_divbfix(ib)) cycle
5018  if(any(typeboundary(:,ib)==bc_special)) then
5019  ! MF nonlinear force-free B field extrapolation and data driven
5020  ! require normal B of the first ghost cell layer to be untouched by
5021  ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
5022  select case (idims)
5023  {case (^d)
5024  if (iside==2) then
5025  ! maximal boundary
5026  ixomin^dd=ixghi^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixglo^dd;
5027  ixomax^dd=ixghi^dd;
5028  else
5029  ! minimal boundary
5030  ixomin^dd=ixglo^dd;
5031  ixomax^dd=ixglo^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixghi^dd;
5032  end if \}
5033  end select
5034  call fixdivb_boundary(ixg^ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
5035  end if
5036  end do
5037  end do
5038 
5039  end subroutine twofl_boundary_adjust
5040 
5041  subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
5043 
5044  integer, intent(in) :: ixG^L,ixO^L,iB
5045  double precision, intent(inout) :: w(ixG^S,1:nw)
5046  double precision, intent(in) :: x(ixG^S,1:ndim)
5047 
5048  double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
5049  integer :: ix^D,ixF^L
5050 
5051  select case(ib)
5052  case(1)
5053  ! 2nd order CD for divB=0 to set normal B component better
5054  {^iftwod
5055  ixfmin1=ixomin1+1
5056  ixfmax1=ixomax1+1
5057  ixfmin2=ixomin2+1
5058  ixfmax2=ixomax2-1
5059  if(slab_uniform) then
5060  dx1x2=dxlevel(1)/dxlevel(2)
5061  do ix1=ixfmax1,ixfmin1,-1
5062  w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
5063  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5064  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5065  enddo
5066  else
5067  do ix1=ixfmax1,ixfmin1,-1
5068  w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
5069  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
5070  +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5071  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5072  -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5073  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5074  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5075  end do
5076  end if
5077  }
5078  {^ifthreed
5079  ixfmin1=ixomin1+1
5080  ixfmax1=ixomax1+1
5081  ixfmin2=ixomin2+1
5082  ixfmax2=ixomax2-1
5083  ixfmin3=ixomin3+1
5084  ixfmax3=ixomax3-1
5085  if(slab_uniform) then
5086  dx1x2=dxlevel(1)/dxlevel(2)
5087  dx1x3=dxlevel(1)/dxlevel(3)
5088  do ix1=ixfmax1,ixfmin1,-1
5089  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5090  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5091  +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5092  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5093  +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5094  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5095  end do
5096  else
5097  do ix1=ixfmax1,ixfmin1,-1
5098  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5099  ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5100  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5101  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5102  +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5103  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5104  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5105  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5106  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5107  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5108  +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5109  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5110  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5111  -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5112  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5113  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5114  /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5115  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5116  end do
5117  end if
5118  }
5119  case(2)
5120  {^iftwod
5121  ixfmin1=ixomin1-1
5122  ixfmax1=ixomax1-1
5123  ixfmin2=ixomin2+1
5124  ixfmax2=ixomax2-1
5125  if(slab_uniform) then
5126  dx1x2=dxlevel(1)/dxlevel(2)
5127  do ix1=ixfmin1,ixfmax1
5128  w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
5129  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
5130  w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
5131  enddo
5132  else
5133  do ix1=ixfmin1,ixfmax1
5134  w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
5135  w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
5136  -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
5137  block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
5138  +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
5139  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
5140  /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
5141  end do
5142  end if
5143  }
5144  {^ifthreed
5145  ixfmin1=ixomin1-1
5146  ixfmax1=ixomax1-1
5147  ixfmin2=ixomin2+1
5148  ixfmax2=ixomax2-1
5149  ixfmin3=ixomin3+1
5150  ixfmax3=ixomax3-1
5151  if(slab_uniform) then
5152  dx1x2=dxlevel(1)/dxlevel(2)
5153  dx1x3=dxlevel(1)/dxlevel(3)
5154  do ix1=ixfmin1,ixfmax1
5155  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5156  w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
5157  -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
5158  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
5159  -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
5160  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
5161  end do
5162  else
5163  do ix1=ixfmin1,ixfmax1
5164  w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
5165  ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
5166  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
5167  block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
5168  -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
5169  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
5170  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
5171  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
5172  w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
5173  block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
5174  -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
5175  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
5176  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
5177  +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
5178  w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
5179  block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
5180  /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
5181  w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
5182  end do
5183  end if
5184  }
5185  case(3)
5186  {^iftwod
5187  ixfmin1=ixomin1+1
5188  ixfmax1=ixomax1-1
5189  ixfmin2=ixomin2+1
5190  ixfmax2=ixomax2+1
5191  if(slab_uniform) then
5192  dx2x1=