DISPATCH
non_ideal_mod.f90
1 !===============================================================================
2 !> non-ideal MHD module
3 !===============================================================================
5  USE const
6  USE io_mod
7  USE io_unit_mod
8  USE link_mod
9  use hydro_parameters
10  USE trace_mod
11  USE mpi_mod
12  USE omp_mod
13  USE index_mod
14  USE mesh_mod
15  implicit none
16  private
17  type, public :: non_ideal_t
18  real:: eta_ohm, etamd, gamma_ad, rho_ions
19  real:: rho0
20  real:: crho
21  real:: t0
22  real:: ct
23  real:: b0
24  real:: cb
25  real:: e0
26  real:: ce
27  real:: mu
28  real:: mp
29  real:: kb
30  real:: na
31  real,dimension(200) :: rhoarr
32  real,dimension(60) :: tarr
33  real,dimension(150) :: barr
34  real,dimension(60) :: earr
35  real,dimension(200,60,150) :: eta_ohm_tbl, eta_ad_tbl
36  real,pointer,dimension(:,:,:,:,:) :: flux
37  logical:: mhd_ad, mhd_ohm, ntestdadm
38  logical:: is_used
39  contains
40  procedure:: init
41  procedure:: non_ideal_comp
42  procedure:: flux_upd
43  procedure:: non_ideal_emf_up
44  procedure:: computejb2
45  procedure:: computambip
46  procedure:: computdifmag
47  procedure:: eta_table
48  procedure:: eta_table_init
49 ! procedure:: comp_weight
50 ! procedure:: read_table
51  end type
52  type(non_ideal_t), public:: non_ideal
53 CONTAINS
54 
55 !===============================================================================
56 !> Read parameters
57 !===============================================================================
58 SUBROUTINE init (self)
59  class(non_ideal_t):: self
60  integer:: iostat
61  real:: gamma_ad=75., rho_ions=1.0, eta_ohm=0.0, etamd=0.0
62  logical:: mhd_ad=.false., mhd_ohm=.false., ntestdadm=.true.
63  logical, save:: first_time=.true.
64  namelist /non_ideal_params/ gamma_ad, eta_ohm, etamd, rho_ions, mhd_ad, mhd_ohm, ntestdadm
65  !---------------------------------------------------------------------------
66  !$omp critical (init_cr)
67  if (first_time) then
68  first_time = .false.
69  rewind(io_unit%input)
70  read (io_unit%input, non_ideal_params,iostat=iostat)
71  if (io%master) write (*, non_ideal_params)
72  end if
73  !$omp end critical (init_cr)
74  self%mhd_Ohm = mhd_ohm
75  self%eta_Ohm = eta_ohm
76  self%mhd_AD = mhd_ad
77  self%gamma_AD = gamma_ad
78  self%rho_ions = rho_ions
79  self%etaMD = etamd
80  self%is_used = mhd_ohm .or. mhd_ad
81  self%ntestDADM=ntestdadm
82  if ( (self%is_used .eqv. .true.) .and. (self%ntestDADM .eqv. .false.) ) then
83  self%eta_Ohm_tbl = 0.d0
84  self%eta_AD_tbl = 0.d0
85  self%kB = 1.3807e-16 ! Boltzmann constant
86  self%mu = 2.37 ! mean molecular weight
87  self%mp = 1.6737236e-24 ! Proton mass
88  self%NA = 6.022140857e23 ! Avogadro constant
89  call eta_table_init(self,self%eta_Ohm_tbl,self%eta_AD_tbl)
90  endif
91 
92 END SUBROUTINE init
93 
94 subroutine non_ideal_comp(self,mesh,idx,bf,uin,qin,fluxambdiff,emfambdiff,u_max,ds,dt,gamma)
95  class(non_ideal_t):: self
96  class(mesh_t):: mesh(:)
97  class(index_t):: idx
98  integer:: n(4), l(3), u(3)
99  real :: u_max
100  real(8) :: gamma
101  real(8) :: ds(3), dt
102  real, dimension(:,:,:,:):: bf, uin, qin, fluxambdiff, emfambdiff
103  ! Output electromotive force
104  real(8), allocatable, dimension(:,:,:)::emfx
105  real(8), allocatable, dimension(:,:,:)::emfy
106  real(8), allocatable, dimension(:,:,:)::emfz
107 
108  ! Output courant vector in the cell
109  real(8), allocatable, dimension(:,:,:,:,:)::bmagij
110  real(8), allocatable, dimension(:,:,:)::jcentersquare,jxbsquare
111  real(8), allocatable, dimension(:,:,:,:)::bemfx,bemfy,bemfz
112  real(8), allocatable, dimension(:,:,:,:)::jemfx,jemfy,jemfz
113  real(8), allocatable, dimension(:,:,:,:)::florentzx,florentzy,florentzz
114  real(8), allocatable, dimension(:,:,:,:)::fluxmd,fluxh,fluxad
115  real(8), allocatable, dimension(:,:,:,:)::emfohmdiss,fluxohm
116 
117  n = shape(uin)
118  allocate (emfx(n(1),n(2),n(3)))
119  allocate (emfy(n(1),n(2),n(3)))
120  allocate (emfz(n(1),n(2),n(3)))
121  allocate (bmagij(n(1),n(2),n(3),3,3))
122  allocate (jcentersquare(n(1),n(2),n(3)))
123  allocate (jxbsquare(n(1),n(2),n(3)))
124  allocate (bemfx(n(1),n(2),n(3),3))
125  allocate (bemfy(n(1),n(2),n(3),3))
126  allocate (bemfz(n(1),n(2),n(3),3))
127  allocate (jemfx(n(1),n(2),n(3),3))
128  allocate (jemfy(n(1),n(2),n(3),3))
129  allocate (jemfz(n(1),n(2),n(3),3))
130  allocate (florentzx(n(1),n(2),n(3),3))
131  allocate (florentzy(n(1),n(2),n(3),3))
132  allocate (florentzz(n(1),n(2),n(3),3))
133  allocate (fluxmd(n(1),n(2),n(3),3))
134  allocate (fluxh(n(1),n(2),n(3),3))
135  allocate (fluxad(n(1),n(2),n(3),3))
136  allocate (emfohmdiss(n(1),n(2),n(3),3))
137  allocate (fluxohm(n(1),n(2),n(3),3))
138 
139  emfx=0.d0
140  emfy=0.d0
141  emfz=0.d0
142  bmagij=0.d0
143  jcentersquare=0.d0
144  jxbsquare=0d0
145  bemfx=0.d0
146  bemfy=0.d0
147  bemfz=0.d0
148  jemfx=0.d0
149  jemfy=0.d0
150  jemfz=0.d0
151  florentzx=0.d0
152  florentzy=0.d0
153  florentzz=0.d0
154  fluxmd=0.d0
155  fluxh=0.d0
156  fluxad=0.d0
157  emfohmdiss=0.d0
158  fluxohm=0.d0
159 
160 
161  if((self%mhd_AD).or.(self%mhd_Ohm)) then
162  ! compute Lorentz Force with current
163  call computejb2(self,bf,qin,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,ds)
164  endif
165 
166 
167  ! AMBIPOLAR DIFFUSION
168 
169  if(self%mhd_AD) then
170  call computambip(self,uin,qin,bemfx,bemfy,bemfz,florentzx,florentzy,florentzz,fluxad,bmagij,emfambdiff,fluxambdiff,jxbsquare,ds,dt,gamma)
171  l = mesh%lo
172  u = mesh%uo
173  u_max = max(u_max, maxval( &
174  uin(l(1):u(1),l(2):u(2),l(3):u(3),idx%bx)**2 + &
175  uin(l(1):u(1),l(2):u(2),l(3):u(3),idx%by)**2 + &
176  uin(l(1):u(1),l(2):u(2),l(3):u(3),idx%bz)**2 / &
177  (self%gamma_AD * self%rho_ions * uin(:,:,:,1) * minval(ds)) ) )
178  endif
179 
180 
181  ! OHMIC DISSIPATION
182 
183  if(self%mhd_Ohm) then
184  call computdifmag(self,uin,qin,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,fluxmd,emfohmdiss,fluxohm,jcentersquare,ds,dt)
185  u_max = max(u_max,self%eta_Ohm/ minval(ds))
186  endif
187 
188 
189  deallocate(emfx,emfy,emfz,bmagij,jcentersquare,jxbsquare,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,emfohmdiss,fluxohm)
190 
191 end subroutine
192 
193 subroutine flux_upd(self,flux,fluxambdiff, dx, dt, idim)
194  class(non_ideal_t):: self
195 
196  integer:: n(5)
197  integer :: i, j, k, idim
198  real, dimension(:,:,:,:,:):: flux
199  real, dimension(:,:,:,:)::fluxambdiff
200 
201  real, allocatable, dimension(:,:,:,:)::fluxohm
202  real(8) :: dx, dt
203 
204 
205  n = shape(flux)
206  allocate (fluxohm(n(1),n(2),n(3),3))
207 
208  fluxohm=0d0
209  if (self%mhd_AD .eqv. .true. .or. self%mhd_Ohm .eqv. .true.) then
210  flux(:,:,:,5,idim) = flux(:,:,:,5,idim) + (fluxambdiff(:,:,:,idim)+fluxohm(:,:,:,idim) ) * dt / dx
211  end if
212 
213  deallocate(fluxohm)
214 
215 end subroutine
216 
217 
218 subroutine non_ideal_emf_up(self,emfup,emfambdiff, dx, dt, idim)
219  class(non_ideal_t):: self
220  integer:: n(4)
221  integer :: idim
222  real, dimension(:,:,:)::emfup
223  real, dimension(:,:,:,:)::emfambdiff
224  real, allocatable, dimension(:,:,:,:)::emfohmdiss
225  real(8) :: dx, dt
226 
227  n = shape(emfambdiff)
228  allocate (emfohmdiss(n(1),n(2),n(3),n(4)))
229  emfohmdiss = 0d0
230 
231  if (self%mhd_AD .eqv. .true. .or. self%mhd_Ohm .eqv. .true.) then
232 
233 #if NDIM==1
234 
235  emfup(:,:,:) = ( emfambdiff(:,:,:,idim) + emfohmdiss(:,:,:,idim) ) * dt / dx
236 #else
237  emfup(2:n(1)-1,2:n(2)-1,2:n(3)-1) = emfup(2:n(1)-1,2:n(2)-1,2:n(3)-1) + &
238  ( emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,idim) + emfohmdiss(2:n(1)-1,2:n(2)-1,2:n(3)-1,idim) ) * dt / dx
239 #endif
240  end if
241 
242  deallocate(emfohmdiss)
243 
244 end subroutine
245 
246 
247 subroutine computejb2(self,bf,q,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,ds)
248  class(non_ideal_t):: self
249 
250  integer :: i, j, k, ibegin, iend, jbegin, jend, kbegin, kend
251  real, dimension(:,:,:,:)::bf
252  real, dimension(:,:,:,:)::q
253 
254 
255  real(8), dimension(:,:,:,:)::bemfx,bemfy,bemfz
256  real(8), dimension(:,:,:,:)::jemfx,jemfy,jemfz
257  real(8), dimension(:,:,:,:,:)::bmagij
258  real(8), dimension(:,:,:,:)::florentzx,florentzy,florentzz
259  real(8), dimension(:,:,:,:)::fluxmd,fluxh,fluxad
260 ! real,dimension(1:self%gn(1),1:self%gn(2),1:self%gn(3),1:3)::jcell
261 
262 
263  ! declare local variables
264  integer :: m
265 
266  integer:: n(4)
267  real(8), allocatable, dimension(:,:,:,:)::bmagijbis
268  real(8), allocatable, dimension(:,:,:,:,:)::jface
269  real(8), allocatable, dimension(:,:,:,:)::bcenter
270  real(8), allocatable, dimension(:,:,:,:,:)::fluxbis,fluxter,fluxquat
271  real(8)::ds(3),dx,dy,dz
272 
273  n = shape(bemfx)
274  allocate (bmagijbis(n(1),n(2),n(3),n(4)))
275  allocate (jface(n(1),n(2),n(3),n(4),3))
276  allocate (bcenter(n(1),n(2),n(3),n(4)))
277  allocate (fluxbis(n(1),n(2),n(3),n(4),3))
278  allocate (fluxter(n(1),n(2),n(3),n(4),3))
279  allocate (fluxquat(n(1),n(2),n(3),n(4),3))
280 
281  bmagijbis = 0.d0
282  jface = 0.d0
283  bcenter = 0.d0
284  fluxbis = 0.d0
285  fluxter = 0.d0
286  fluxquat = 0.d0
287 
288  dx = ds(1) ; dy = ds(2) ; dz = ds(3)
289  ! Compute before the loop because bmagijbis has to be computed with index i+1/j+1/k+1 in loop over i/j/k
290  bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)= 0.25d0* (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) + bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) + bf(2:n(1)-1,1:n(2)-2,1:n(3)-2,1) )
291  bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)= 0.25d0* (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) + bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,2) + bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,2) + bf(1:n(1)-2,2:n(2)-1,1:n(3)-2,2) )
292  bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)= 0.25d0* (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) + bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,3) + bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,3) + bf(1:n(1)-2,1:n(2)-2,2:n(3)-1,3) )
293  ! Define the 1 and 2 components of bcenter, since they are needed further below, when computing the current on the cell interfaces (jface)
294  ! magnetic field at center of cells
295  bcenter(1:n(1),1:n(2),1:n(3),1:3)=q(1:n(1),1:n(2),1:n(3),6:8)
296 
297  !!!!!!!!!!!!!!!!!!
298  ! EMF x
299  !!!!!!!!!!!!!!!!!!
300  bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=0.25d0* ( q(2:n(1)-1,2:n(2)-1,2:n(3)-1,6) + q(2:n(1)-1,1:n(2)-2,2:n(3)-1,6) + q(2:n(1)-1,2:n(2)-1,1:n(3)-2,6) + q(2:n(1)-1,1:n(2)-2,1:n(3)-2,6) )
301  bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=0.5d0 * ( bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) + bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,2) )
302  bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=0.5d0 * ( bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) + bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,3) )
303  !!!!!!!!!!!!!!!!!!
304  ! EMF y
305  !!!!!!!!!!!!!!!!!!
306  bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=0.5d0 * ( bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) )
307  bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=0.25d0* ( q(2:n(1)-1,2:n(2)-1,2:n(3)-1,7) + q(1:n(1)-2,2:n(2)-1,2:n(3)-1,7) + q(2:n(1)-1,2:n(2)-1,1:n(3)-2,7) + q(1:n(1)-2,2:n(2)-1,1:n(3)-2,7) )
308  bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=0.5d0 * ( bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,3) + bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) )
309  !!!!!!!!!!!!!!!!!!
310  ! EMF z
311  !!!!!!!!!!!!!!!!!!
312  bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=0.5d0 * ( bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) )
313  bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=0.5d0 * ( bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) + bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,2) )
314  bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=0.25d0* ( q(2:n(1)-1,2:n(2)-1,2:n(3)-1,8) + q(1:n(1)-2,2:n(2)-1,2:n(3)-1,8) + q(2:n(1)-1,1:n(2)-2,2:n(3)-1,8) + q(1:n(1)-2,1:n(2)-2,2:n(3)-1,8) )
315  do m=1,3
316  ! bmagij is the value of the magnetic field Bi where Bj
317  ! is naturally defined; Ex bmagij(i,j,k,1,2) is Bx at i,j-1/2,k
318  ! and we can write it Bx,y
319  !! m+5 mandatory cf Bx=uin(i,j,k,6)
320  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,m,m)=bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,m)
321  end do
322  ! case Bx,y
323  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,2) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,6) + q(2:n(1)-1,1:n(2)-2,2:n(3)-1,6) )
324  ! case Bx,z
325  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,6) + q(2:n(1)-1,2:n(2)-1,1:n(3)-2,6) )
326  ! case By,y
327  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,7) + q(1:n(1)-2,2:n(2)-1,2:n(3)-1,7) )
328  ! case By,z
329  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,3) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,7) + q(2:n(1)-1,2:n(2)-1,1:n(3)-2,7) )
330  ! case Bz,x
331  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,8) + q(1:n(1)-2,2:n(2)-1,2:n(3)-1,8) )
332  ! case Bz,y
333  bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2) = 0.5d0 * (q(2:n(1)-1,2:n(2)-1,2:n(3)-1,8) + q(2:n(1)-1,1:n(2)-2,2:n(3)-1,8) )
334  jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) = (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,3) ) / dy - (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,2) ) / dz
335  jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) = (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,2) - bmagij(2:n(1)-1,2:n(2)-1,1:n(3)-2,1,2) ) / dz - (bmagijbis(1:n(1)-2,2:n(2)-1,2:n(3)-1,3) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dx
336  jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) = (bmagijbis(1:n(1)-2,2:n(2)-1,2:n(3)-1,2) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) ) / dx - (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3)- bmagij(2:n(1)-1,1:n(2)-2,2:n(3)-1,1,3) ) / dy
337  jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) = (bmagijbis(2:n(1)-1,3:n(2) ,2:n(3)-1,3) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dy - (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1)- bmagij(2:n(1)-1,2:n(2)-1,1:n(3)-2,2,1) ) / dz
338  jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) = (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - bf(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) ) / dz - (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,3) ) / dx
339  jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) = (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,3) - bmagij(1:n(1)-2,2:n(2)-1,2:n(3)-1,2,3) ) / dx - (bmagijbis(2:n(1)-1,3:n(2) ,2:n(3)-1,1) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dy
340  jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) = (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1) - bmagij(2:n(1)-1,1:n(2)-2,2:n(3)-1,3,1) ) / dy - (bmagijbis(2:n(1)-1,2:n(2)-1,3:n(3) ,2) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) ) / dz
341  jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) = (bmagijbis(2:n(1)-1,2:n(2)-1,3:n(3) ,1) - bmagijbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dz - (bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2) - bmagij(1:n(1)-2,2:n(2)-1,2:n(3)-1,3,2) ) / dx
342  jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) = (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - bf(1:n(1)-2,2:n(2)-1,2:n(3)-1,2) ) / dx - (bf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - bf(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) ) / dy
343 
344 
345 if(self%mhd_AD) then
346 
347 ! EMF x
348  florentzx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
349  florentzx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
350  florentzx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
351  florentzy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
352  florentzy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
353  florentzy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
354  florentzz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
355  florentzz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
356  florentzz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
357 
358 
359 endif
360 
361 
362 ! computation of current on faces
363 
364 if((self%mhd_AD).or.(self%mhd_Ohm)) then
365 
366 ! face at i-1/2,j,k
367 
368  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,1)= (bemfz(2:n(1)-1,3:n(2) ,2:n(3)-1,3) - bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dy - (bemfy(2:n(1)-1,2:n(2)-1,3:n(3) ,2) - bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) ) / dz
369  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1)= (bemfy(2:n(1)-1,2:n(2)-1,3:n(3) ,1) - bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dz - (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - bcenter(1:n(1)-2,2:n(2)-1,2:n(3)-1,3) ) / dx
370  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1)= (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - bcenter(1:n(1)-2,2:n(2)-1,2:n(3)-1,2) ) / dx - (bemfz(2:n(1)-1,3:n(2) ,2:n(3)-1,1) - bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dy
371  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,2)= (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - bcenter(2:n(1)-1,2:n(2)-1,1:n(3)-2,3) ) / dy - (bemfx(2:n(1)-1,2:n(2)-1,3:n(3) ,2) - bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) ) / dz
372  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,2)= (bemfx(2:n(1)-1,2:n(2)-1,3:n(3) ,1) - bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dz - (bemfz(3:n(1) ,2:n(2)-1,2:n(3)-1,3) - bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dx
373  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2)= (bemfz(3:n(1) ,2:n(2)-1,2:n(3)-1,2) - bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dx - (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - bcenter(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) ) / dy
374  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3)= (bemfx(2:n(1)-1,3:n(2) ,2:n(3)-1,3) - bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dy - (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - bcenter(2:n(1)-1,2:n(2)-1,1:n(3)-2,2) ) / dz
375  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,3)= (bcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - bcenter(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) ) / dz - (bemfy(3:n(1) ,2:n(2)-1,2:n(3)-1,3) - bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) ) / dx
376  jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,3)= (bemfy(3:n(1) ,2:n(2)-1,2:n(3)-1,2) - bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) ) / dx - (bemfx(2:n(1)-1,3:n(2) ,2:n(3)-1,1) - bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) ) / dy
377 
378  fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:) = jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:) - jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)
379  fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:) = jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:) - jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)
380  fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:) = jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:) - jface(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)
381 
382  fluxmd(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,1)
383  fluxmd(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,2)
384  fluxmd(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,3)
385 
386 
387 endif
388 
389 if(self%mhd_AD) then
390 
391  fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:) - fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)
392  fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:) - fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)
393  fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)=fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:) - fluxbis(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)
394 
395  fluxh(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,1)
396  fluxh(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,2)
397  fluxh(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,3)
398 
399  fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:) - fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)
400  fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:) - fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)
401  fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,:)=fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:) - fluxter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,:)*bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,:)
402 
403  fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,1)
404  fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,2)
405  fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=fluxquat(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,3)
406 
407 
408 endif
409 
410  deallocate (bmagijbis, jface, bcenter, fluxbis, fluxter, fluxquat)
411 
412 
413 
414 end subroutine computejb2
415 !###########################################################
416 !###########################################################
417 !###########################################################
418 !###########################################################
419 subroutine computdifmag(self,u,q,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,fluxmd,emfohmdiss,fluxohm,jcentersquare,ds,dt)
420  class(non_ideal_t):: self
421 
422  integer :: i, j, k
423 
424  integer:: n(4)
425  real, dimension(:,:,:,:)::u
426  real, dimension(:,:,:,:)::q
427 
428  real(8), dimension(:,:,:,:)::bemfx,bemfy,bemfz
429  real(8), dimension(:,:,:,:)::jemfx,jemfy,jemfz
430  real(8), dimension(:,:,:,:,:)::bmagij
431  real(8), dimension(:,:,:,:)::fluxmd
432 
433  real(8), dimension(:,:,:,:):: emfohmdiss,fluxohm
434  real(8), dimension(:,:,:)::jcentersquare
435 
436  ! declare local variables
437  integer :: m, h
438 
439 
440  ! WARNING following quantities defined with three components even
441  ! if ndim<3 !
442  real(8), allocatable, dimension(:,:,:,:)::jcenter
443  real(8), allocatable, dimension(:,:,:,:,:)::jface
444  real(8), allocatable, dimension(:,:,:,:)::jemf
445 
446  real(8)::rhocell,bcell,cv
447 
448  real(8), allocatable, dimension(:,:,:)::bsquarex,bsquarey,bsquarez,dtlim
449  real(8), allocatable, dimension(:,:,:)::rhox,rhoy,rhoz,epsx,epsy,epsz
450  real(8), allocatable, dimension(:,:,:)::rhof,pf,bsqf,epsf
451 
452  integer, dimension(3) :: index_i,index_j,index_k
453 
454  real(8):: ds(3)
455  real(8):: dx, dy, dz, dt
456 
457  n = shape(u)
458  allocate (jcenter(n(1),n(2),n(3),3))
459  allocate (jface(n(1),n(2),n(3),3,3))
460  allocate (jemf(n(1),n(2),n(3),3))
461  allocate (bsquarex(n(1)-2,n(2)-2,n(3)-2) )
462  allocate (bsquarey(n(1)-2,n(2)-2,n(3)-2) )
463  allocate (bsquarez(n(1)-2,n(2)-2,n(3)-2) )
464  allocate (dtlim(n(1)-2,n(2)-2,n(3)-2) )
465  allocate (rhox(n(1)-2,n(2)-2,n(3)-2) )
466  allocate (rhoy(n(1)-2,n(2)-2,n(3)-2) )
467  allocate (rhoz(n(1)-2,n(2)-2,n(3)-2) )
468  allocate (epsx(n(1)-2,n(2)-2,n(3)-2) )
469  allocate (epsy(n(1)-2,n(2)-2,n(3)-2) )
470  allocate (epsz(n(1)-2,n(2)-2,n(3)-2) )
471  allocate (rhof(n(1)-2,n(2)-2,n(3)-2) )
472  allocate (pf(n(1)-2,n(2)-2,n(3)-2) )
473  allocate (bsqf(n(1)-1,n(2)-1,n(3)-1) )
474  allocate (epsf(n(1)-1,n(2)-1,n(3)-1) )
475 
476  jcenter = 0.d0
477  jface = 0.d0
478  jemf = 0.d0
479  bsquarex = 0.d0
480  bsquarey = 0.d0
481  bsquarez = 0.d0
482  dtlim = 0.d0
483  rhox = 0.d0
484  rhoy = 0.d0
485  rhoz = 0.d0
486  epsx = 0.d0
487  epsy = 0.d0
488  epsz = 0.d0
489  rhof = 0.d0
490  pf = 0.d0
491  bsqf = 0.d0
492  epsf = 0.d0
493 
494  index_i = (/1,0,0/)
495  index_j = (/0,1,0/)
496  index_k = (/0,0,1/)
497 
498  dx = ds(1); dy = ds(2) ; dz = ds(3)
499  dtlim = dt !neil
500 
501  jemf(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=jemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
502  jemf(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=jemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
503  jemf(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=jemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
504  rhox=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) + u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) + u(2:n(1)-1,1:n(2)-2,1:n(3)-2,1))
505  rhoy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) + u(1:n(1)-2,2:n(2)-1,1:n(3)-2,1))
506  rhoz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) + u(1:n(1)-2,1:n(2)-2,2:n(3)-1,1))
507  !epsx=0.25d0*(u(2:ng-1,2:ng-1,2:ng-1,nvar) + u(2:ng-1,j-1,2:ng-1,nvar) + u(2:ng-1,2:ng-1,1:ng-2,nvar) + u(2:ng-1,1:ng-2,1:ng-2,nvar))
508  !epsy=0.25d0*(u(2:ng-1,2:ng-1,2:ng-1,nvar) + u(i-1,2:ng-1,2:ng-1,nvar) + u(2:ng-1,2:ng-1,1:ng-2,nvar) + u(1:ng-2,2:ng-1,1:ng-2,nvar))
509  !epsz=0.25d0*(u(2:ng-1,2:ng-1,2:ng-1,nvar) + u(i-1,2:ng-1,2:ng-1,nvar) + u(2:ng-1,1:ng-2,2:ng-1,nvar) + u(1:ng-2,1:ng-2,2:ng-1,nvar))
510  if(self%mhd_Ohm)then
511  bsquarex=bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
512  bsquarey=bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
513  bsquarez=bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
514  end if
515 
516  ! WARNING dB/dt=-curl(eta*J)
517  emfohmdiss(2:n(1)-1,2:n(2)-1,2:n(3)-1,:)=-self%etaMD*jemf(2:n(1)-1,2:n(2)-1,2:n(3)-1,:)
518  !!!!!!!!!!!!!!!!!!!!!!!
519  !
520  ! compute j at center of cells
521  !
522  ! mandatory for non isotherm case
523 
524  ! bmagij is the value of the magnetic field Bi where Bj
525  ! is naturally defined; Ex bmagij(i,j,k,1,2) is Bx at i,j-1/2,k
526  ! and we can write it Bx,y
527 
528  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)= (bmagij(2:n(1)-1,3:n(2) ,2:n(3)-1,3,2) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2) ) / dz - (bmagij(2:n(1)-1,2:n(2)-1,3:n(3) ,2,3) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,3) ) / dy
529  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)= (bmagij(2:n(1)-1,2:n(2)-1,3:n(3) ,1,3) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3) ) / dx - (bmagij(3:n(1) ,2:n(2)-1,2:n(3)-1,3,1) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1) ) / dz
530  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)= (bmagij(3:n(1) ,2:n(2)-1,2:n(3)-1,2,1) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1) ) / dy - (bmagij(2:n(1)-1,3:n(2) ,2:n(3)-1,2,1) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1) ) / dx
531 
532  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
533  !
534  do h = 1,3
535 
536  rhof=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(2-index_i(h):n(1)-index_i(h),2-index_j(h):n(2)-index_j(h),2-index_k(h):n(3)-index_k(h),1))
537  !epsf=0.5d0*(u(2:ng-1,2:ng-1,2:ng-1,nvar)+u(2-index_i(h):ng-index_i(h),2-index_j(h):ng-index_j(h),2-index_k(h):ng-index_k(h),nvar))
538  bsqf=bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,h)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,h)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,h)**2
539 
540  fluxohm(2:n(1)-1,2:n(2)-1,2:n(3)-1,h)=self%etaMD*fluxmd(2:n(1)-1,2:n(2)-1,2:n(3)-1,h)
541 
542  enddo
543 
544  deallocate(jcenter, jface, jemf, bsquarex, bsquarey, bsquarez, dtlim, rhox, rhoy, rhoz, epsx, epsy, epsz, rhof, pf, bsqf, epsf )
545 
546 end subroutine computdifmag
547 !###########################################################
548 !###########################################################
549 !###########################################################
550 !###########################################################
551 SUBROUTINE computambip(self,u,q,bemfx,bemfy,bemfz,florentzx,florentzy,florentzz,fluxad,bmagij,emfambdiff,fluxambdiff,jxbsquare,ds,dt,gamma)
552  class(non_ideal_t):: self
553 
554  real, dimension(:,:,:,:)::u
555  real, dimension(:,:,:,:)::q
556 
557  real(8), dimension(:,:,:,:)::bemfx,bemfy,bemfz
558  real(8), dimension(:,:,:,:)::florentzx,florentzy,florentzz
559  real(8), dimension(:,:,:,:)::fluxad
560  real(8), dimension(:,:,:,:,:)::bmagij
561  real, dimension(:,:,:,:)::emfambdiff
562  real, dimension(:,:,:,:)::fluxambdiff
563  real(8), dimension(:,:,:)::jxbsquare
564 
565 ! declare local variables
566  integer :: m, ivar
567  integer:: n(4)
568 
569  real(8), allocatable, dimension(:,:,:)::rhocellmin
570  real(8), allocatable, dimension(:,:,:)::rhofx,rhofy,rhofz
571  real(8), allocatable, dimension(:,:,:)::bsquarex,bsquarey,bsquarez,bsquare
572  real(8), allocatable, dimension(:,:,:)::bsquarexx,bsquareyy,bsquarezz
573  real(8), allocatable, dimension(:,:,:)::betaad2,betaad,eta_ad,eta_ohm
574  real(8), allocatable, dimension(:,:,:)::rhox,rhoy,rhoz,rhocell,bcell,bcellold
575  real(8)::dtlim,cv
576 
577  real(8), allocatable, dimension(:,:,:,:)::florentz
578  real(8), allocatable, dimension(:,:,:) ::bsquaremax
579 
580  real(8), allocatable, dimension(:,:,:,:)::jcenter
581  real(8), allocatable, dimension(:,:,:,:)::jxb
582 
583  real(8), allocatable, dimension(:,:,:)::uxfx,uxfy,uxfz,uyfx,uyfy,uyfz,uzfx,uzfy,uzfz,uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,brms2,eth,brms,rho,b,p,t
584  real(8), allocatable, dimension(:,:,:)::etfx,etfy,etfz,etx,ety,etz
585 
586  real(8) :: ds(3)
587  real(8) :: dx, dy, dz, dt, c2_4pi, mp, kb, bcode_to_gauss, pi, scale_v, scale_d, scale_l, scale_t, scale_t2
588  real(8) :: gamma
589 
590  n = shape(u)
591 
592  allocate (rhocellmin(n(1),n(2),n(3)))
593  allocate (rhofx(n(1)-2,n(2)-2,n(3)-2))
594  allocate (rhofy(n(1)-2,n(2)-2,n(3)-2))
595  allocate (rhofz(n(1)-2,n(2)-2,n(3)-2))
596  allocate (bsquarex(n(1)-2,n(2)-2,n(3)-2))
597  allocate (bsquarey(n(1)-2,n(2)-2,n(3)-2))
598  allocate (bsquarez(n(1)-2,n(2)-2,n(3)-2))
599  allocate (bsquarexx(n(1)-2,n(2)-2,n(3)-2))
600  allocate (bsquareyy(n(1)-2,n(2)-2,n(3)-2))
601  allocate (bsquarezz(n(1)-2,n(2)-2,n(3)-2))
602  allocate (betaad(n(1)-2,n(2)-2,n(3)-2))
603  allocate (betaad2(n(1)-2,n(2)-2,n(3)-2))
604  allocate (eta_ad(n(1)-2,n(2)-2,n(3)-2))
605  allocate (eta_ohm(n(1)-2,n(2)-2,n(3)-2))
606  allocate (rhox(n(1)-2,n(2)-2,n(3)-2))
607  allocate (rhoy(n(1)-2,n(2)-2,n(3)-2))
608  allocate (rhoz(n(1)-2,n(2)-2,n(3)-2))
609  allocate (rhocell(n(1)-2,n(2)-2,n(3)-2))
610  allocate (bcell(n(1)-2,n(2)-2,n(3)-2))
611  allocate (bcellold(n(1)-2,n(2)-2,n(3)-2))
612  allocate (florentz(n(1),n(2),n(3),3))
613  allocate (bsquaremax(n(1),n(2),n(3)))
614  allocate (jcenter(n(1),n(2),n(3),3))
615  allocate (jxb(n(1),n(2),n(3),3))
616  allocate (uxfx(n(1)-2,n(2)-2,n(3)-2))
617  allocate (uxfy(n(1)-2,n(2)-2,n(3)-2))
618  allocate (uxfz(n(1)-2,n(2)-2,n(3)-2))
619  allocate (uyfx(n(1)-2,n(2)-2,n(3)-2))
620  allocate (uyfy(n(1)-2,n(2)-2,n(3)-2))
621  allocate (uyfz(n(1)-2,n(2)-2,n(3)-2))
622  allocate (uzfx(n(1)-2,n(2)-2,n(3)-2))
623  allocate (uzfy(n(1)-2,n(2)-2,n(3)-2))
624  allocate (uzfz(n(1)-2,n(2)-2,n(3)-2))
625  allocate (uxx(n(1)-2,n(2)-2,n(3)-2))
626  allocate (uxy(n(1)-2,n(2)-2,n(3)-2))
627  allocate (uxz(n(1)-2,n(2)-2,n(3)-2))
628  allocate (uyx(n(1)-2,n(2)-2,n(3)-2))
629  allocate (uyy(n(1)-2,n(2)-2,n(3)-2))
630  allocate (uyz(n(1)-2,n(2)-2,n(3)-2))
631  allocate (uzx(n(1)-2,n(2)-2,n(3)-2))
632  allocate (uzy(n(1)-2,n(2)-2,n(3)-2))
633  allocate (uzz(n(1)-2,n(2)-2,n(3)-2))
634  allocate (etfx(n(1)-2,n(2)-2,n(3)-2))
635  allocate (etfy(n(1)-2,n(2)-2,n(3)-2))
636  allocate (etfz(n(1)-2,n(2)-2,n(3)-2))
637  allocate (etx(n(1)-2,n(2)-2,n(3)-2))
638  allocate (ety(n(1)-2,n(2)-2,n(3)-2))
639  allocate (etz(n(1)-2,n(2)-2,n(3)-2))
640  allocate (brms2(n(1)-2,n(2)-2,n(3)-2))
641  allocate (eth(n(1)-2,n(2)-2,n(3)-2))
642  allocate (brms(n(1)-2,n(2)-2,n(3)-2))
643  allocate (rho(n(1)-2,n(2)-2,n(3)-2))
644  allocate (b(n(1)-2,n(2)-2,n(3)-2))
645  allocate (p(n(1)-2,n(2)-2,n(3)-2))
646  allocate (t(n(1)-2,n(2)-2,n(3)-2))
647 
648  rhocellmin = 0.d0
649  rhofx = 0.d0
650  rhofy = 0.d0
651  rhofz = 0.d0
652  bsquarex = 0.d0
653  bsquarey = 0.d0
654  bsquarez = 0.d0
655  bsquarexx = 0.d0
656  bsquareyy = 0.d0
657  bsquarezz = 0.d0
658  betaad = 0.d0
659  betaad2 = 0.d0
660  eta_ad = 0.d0
661  eta_ohm = 0.d0
662  rhox = 0.d0
663  rhoy = 0.d0
664  rhoz = 0.d0
665  rhocell = 0.d0
666  bcell = 0.d0
667  bcellold = 0.d0
668  florentz = 0.d0
669  bsquaremax = 0.d0
670  jcenter = 0.d0
671  jxb = 0.d0
672 
673  dx = ds(1) ; dy = ds(2) ; dz = ds(3)
674 
675  if ( self%ntestDADM .eqv. .false. ) then
676  scale_v = 1e6
677  scale_d = 1e-17
678  scale_l = 5e15
679  scale_t = scale_l/scale_v
680  scale_t2 = (scale_v)**2 * self%mp / self%kB
681  pi = 4.*atan(1.d0)
682  c2_4pi = 7.152426325648638e+19
683  bcode_to_gauss = (4.d0*pi*scale_d)**0.5*scale_v
684  endif
685 
686  !dt est deja dtnew, qui a été choisi comme le dt normal (avec la condition de courant) ou le dt normal seuillé si le dtAD est trop faible(bricolo)
687  dtlim=dt!*coefalfven
688 
689  jxb=0.0d0
690  jxbsquare=0.0d0
691  jcenter=0.0d0
692 
693  rhox=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) + u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) + u(2:n(1)-1,1:n(2)-2,1:n(3)-2,1))
694  rhoy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1) + u(1:n(1)-2,2:n(2)-1,1:n(3)-2,1))
695  rhoz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) + u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1) + u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1) + u(1:n(1)-2,1:n(2)-2,2:n(3)-1,1))
696 
697  rhofx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1))
698  rhofy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1))
699  rhofz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1))
700 
701  rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)=min(rhox,rhoy,rhoz,rhofx,rhofy,rhofz)
702  rhocell = u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
703 
704  bsquarex=bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
705  bsquarey=bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
706  bsquarez=bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2+bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2+bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
707 
708  bsquarexx=bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,1)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1)**2
709  bsquareyy=bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,2)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,2)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2)**2
710  bsquarezz=bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,3)**2+bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,3)**2
711 
712  bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)=max(bsquarex,bsquarey,bsquarez,bsquarexx,bsquareyy,bsquarezz)
713 
714  ! EMF x
715  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)= florentzx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) - florentzx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
716  rhox=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1)+u(2:n(1)-1,1:n(2)-2,1:n(3)-2,1))
717  bcell = bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)
718  bcellold=bcell
719  rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
720 
721  if (self%ntestDADM) then
722  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhox)
723  else
724  ! Call resistivity table dependent on density, thermal energy and magnetic field strength
725  ! Eth = Etot - kinetic energy - magnetic energy
726  uxx=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,2)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,2)+u(2:n(1)-1,1:n(2)-2,1:n(3)-2,2))
727  uyx=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,3)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,3)+u(2:n(1)-1,1:n(2)-2,1:n(3)-2,3))
728  uzx=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,4)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,4)+u(2:n(1)-1,1:n(2)-2,1:n(3)-2,4))
729  etx=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,5)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,5)+u(2:n(1)-1,1:n(2)-2,1:n(3)-2,5))
730  brms2 = bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2 + bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2 + bemfx(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
731  eth = etx - 0.5 / rhox * ( uxx**2 + uyx**2 + uzx**2 ) - &
732  0.5 * brms2
733  brms = brms2**0.5
734  rho = rhox * scale_d
735  p = eth * (gamma-1.)
736  t = p / rhox * scale_t2 * self%mu
737  b = brms * bcode_to_gauss
738  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
739  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
740  ! from seconds to code time
741 
742  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
743  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
744  betaad2 = eta_ad / brms2
745  endif
746 
747  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*betaad2
748 
749  ! EMF y
750  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)= florentzy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) - florentzy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
751  rhoy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1)+u(1:n(1)-2,2:n(2)-1,1:n(3)-2,1))
752  bcell = bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)
753  bcellold=bcell
754  rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
755  ! comparison with hydro+idealMHD
756 
757  if (self%ntestDADM) then
758  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhoy)
759  else
760  ! Call resistivity table dependent on density, thermal energy and magnetic field strength
761  ! Eth = Etot - kinetic energy - magnetic energy
762  uxy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,2)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,2)+u(1:n(1)-2,2:n(2)-1,1:n(3)-2,2))
763  uyy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,3)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,3)+u(1:n(1)-2,2:n(2)-1,1:n(3)-2,3))
764  uzy=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,4)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,4)+u(1:n(1)-2,2:n(2)-1,1:n(3)-2,4))
765  ety=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,5)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,5)+u(1:n(1)-2,2:n(2)-1,1:n(3)-2,5))
766  brms2 = bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2 + bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2 + bemfy(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
767  eth = ety - 0.5 / rhoy * ( uxy**2 + uyy**2 + uzy**2 ) - &
768  0.5 * brms2
769  brms = brms2**0.5
770  rho = rhoy*scale_d
771  p = eth * (gamma-1.)
772  t = p / rhoy * scale_t2 * self%mu
773  b = brms * bcode_to_gauss
774  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
775  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
776  ! from seconds to code time
777  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
778  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
779  betaad2 = eta_ad / brms2
780  endif
781  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*betaad2
782 
783  ! EMF z
784  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)= florentzz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) - florentzz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
785  rhoz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1)+u(1:n(1)-2,1:n(2)-2,2:n(3)-1,1))
786  bcell = bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)
787  bcellold=bcell
788  rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
789 
790  if (self%ntestDADM) then
791  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhoz)
792  else
793  ! Call resistivity table dependent on density, thermal energy and magnetic field strength
794  ! Eth = Etot - kinetic energy - magnetic energy
795  uxz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,2)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,2)+u(1:n(1)-2,1:n(2)-2,2:n(3)-1,2))
796  uyz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,3)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,3)+u(1:n(1)-2,1:n(2)-2,2:n(3)-1,3))
797  uzz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,4)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,4)+u(1:n(1)-2,1:n(2)-2,2:n(3)-1,4))
798  etz=0.25d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,5)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,5)+u(1:n(1)-2,1:n(2)-2,2:n(3)-1,5))
799  brms2 = bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)**2 + bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2 + bemfz(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2
800  eth = etz - 0.5 / rhoz * ( uxz**2 + uyz**2 + uzz**2 ) - &
801  0.5 * brms2
802  brms = brms2**0.5
803  rho = rhoz*scale_d
804  p = eth * (gamma-1.)
805  t = p / rhoz * scale_t2 * self%mu
806  b = brms * bcode_to_gauss
807  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
808  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
809  ! from seconds to code time
810  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
811  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
812  betaad2 = eta_ad / brms2
813  endif
814 
815  emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=emfambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*betaad2
816 
817  ! energy flux on faces
818  bcell = bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)
819  rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
820  rhofx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,1))
821  if (self%ntestDADM) then
822  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhofx)
823  else
824  uxfx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,2))
825  uyfx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,3))
826  uzfx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,4))
827  etfx=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(1:n(1)-2,2:n(2)-1,2:n(3)-1,5))
828  brms2 = bsquarex
829  eth = etfx - 0.5 / rhofx * ( uxfx**2 + uyfx**2 + uzfx**2 ) - &
830  0.5 * brms2
831  brms = brms2**0.5
832  rho = rhofx*scale_d
833  p = eth * (gamma-1.)
834  t = p / rhofx * scale_t2 *self%mu
835  b = brms * bcode_to_gauss
836  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
837  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
838  ! from seconds to code time
839  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
840  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
841  betaad2 = eta_ad / brms2
842  endif
843  fluxambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)=-betaad2*fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
844 
845  rhofy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,1))
846  if (self%ntestDADM) then
847  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhofy)
848  else
849  uxfy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,2))
850  uyfy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,3))
851  uzfy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,4))
852  etfy=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(2:n(1)-1,1:n(2)-2,2:n(3)-1,5))
853  brms2 = bsquarey
854  eth = etfy - 0.5 / rhofy * ( uxfy**2 + uyfy**2 + uzfy**2 ) - &
855  0.5 * brms2
856  brms = brms2**0.5
857  rho = rhofy*scale_d
858  p = eth * (gamma-1.)
859  t = p / rhofy * scale_t2 * self%mu
860  b = brms * bcode_to_gauss
861  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
862  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
863  ! from seconds to code time
864  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
865  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
866  betaad2 = eta_ad / brms2
867  endif
868  fluxambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)=-betaad2*fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)
869 
870  rhofz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,1))
871  if (self%ntestDADM) then
872  betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhofz)
873  else
874  uxfz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,2))
875  uyfz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,3))
876  uzfz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,4))
877  etfz=0.5d0*(u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5)+u(2:n(1)-1,2:n(2)-1,1:n(3)-2,5))
878  brms2 = bsquarez
879  eth = etfz - 0.5 / rhofz * ( uxfz**2 + uyfz**2 + uzfz**2 ) - &
880  0.5 * brms2
881  brms = brms2**0.5
882  rho = rhofz*scale_d
883  p = eth * (gamma-1.)
884  t = p / rhofz * scale_t2 * self%mu
885  b = brms * bcode_to_gauss
886  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
887  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
888  ! from seconds to code time
889  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
890  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
891  betaad2 = eta_ad / brms2
892  endif
893  fluxambdiff(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)=-betaad2*fluxad(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)
894 
895  bcellold=bcell
896 
897 
898  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)= ( bmagij(2:n(1)-1,3:n(2) ,2:n(3)-1,3,2) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2) ) / dz - ( bmagij(2:n(1)-1,2:n(2)-1,3:n(3) ,2,3) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,2) ) / dy
899  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)= ( bmagij(2:n(1)-1,2:n(2)-1,3:n(3) ,1,3) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,3) ) / dx - ( bmagij(3:n(1) ,2:n(2)-1,2:n(3)-1,3,1) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,3,1) ) / dz
900  jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)= ( bmagij(3:n(1) ,2:n(2)-1,2:n(3)-1,2,1) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,2,1) ) / dy - ( bmagij(2:n(1)-1,3:n(2) ,2:n(3)-1,1,2) - bmagij(2:n(1)-1,2:n(2)-1,2:n(3)-1,1,2) ) / dx
901 
902  jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) = jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,8) - jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,7)
903  jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) = jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,6) - jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,8)
904  jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,3) = jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,7) - jcenter(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,6)
905 
906  if (self%ntestDADM) then
907  betaad=1.d0/(self%gamma_AD*self%rho_ions*u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1))
908  else
909  brms2 = u(2:n(1)-1,2:n(2)-1,2:n(3)-1,6)**2 + u(2:n(1)-1,2:n(2)-1,2:n(3)-1,7)**2 + u(2:n(1)-1,2:n(2)-1,2:n(3)-1,8)**2
910  eth = u(2:n(1)-1,2:n(2)-1,2:n(3)-1,5) - 0.5 / u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) * ( u(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)**2 + u(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)**2 + u(2:n(1)-1,2:n(2)-1,2:n(3)-1,4)**2 ) - &
911  0.5 * brms2
912  brms = brms2**0.5
913  rhocell = u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
914  ! Convert variables to units used in the table
915  ! rho in cm^(-3) ; T in K ; B in Gauss
916  rho = rhocell * scale_d
917  p = eth * (gamma-1.)
918  t = p / rhocell * scale_t2 * self%mu
919  b = brms * bcode_to_gauss
920  call eta_table(self,eta_ad,eta_ohm,rho,t,b)
921  ! Convert back to code units. Reminder: eta_Masson = c^2/(4*pi) * eta_Marchand
922  ! from seconds to code time
923  eta_ad = eta_ad * c2_4pi ! account for necessary prefactor
924  eta_ad = eta_ad / scale_v**2 / scale_t ! transfer back to code units
925  betaad = eta_ad / brms2
926  endif
927  jxbsquare(2:n(1)-1,2:n(2)-1,2:n(3)-1)=(jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)*jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,1) +&
928  & jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,2)*jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,2) +&
929  & jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,3)*jxb(2:n(1)-1,2:n(2)-1,2:n(3)-1,3))*&
930  & betaad*dtlim
931 
932  deallocate (rhocellmin, rhofx, rhofy, rhofz, bsquarex, bsquarey, bsquarez, bsquarexx, bsquareyy, bsquarezz, &
933  betaad, betaad2, eta_ad, eta_ohm, rhox, rhoy, rhoz, rhocell, bcell, bcellold, jcenter, jxb, &
934  uxfx, uxfy, uxfz, uyfx, uyfy, uyfz, uzfx, uzfy, uzfz, uxx, uxy, uxz, uyx, uyy, uyz, uzx, uzy, uzz, &
935  etfx, etfy, etfz, etx, ety, etz, brms2, eth, brms, rho, b, p, t )
936 
937 
938 end SUBROUTINE computambip
939 ! fin modif nimhd
940 
941 
942 subroutine eta_table(self,eta_AD,eta_Ohm,rho,T,B)
943  class(non_ideal_t):: self
944  real(8), dimension(:,:,:) :: rho, t, b
945  real(8), dimension(:,:,:) :: eta_ad, eta_ohm
946  real(8) :: rho0, t0, b0
947  real(8) :: crho, ct, cb
948  integer(8) :: rhoidx, tidx, bidx
949  real :: w1(2), w2(2),w3(2)
950  integer :: n(3),i,j,k,iw,jw,kw,irho,it,ib
951 
952  n =shape(rho)
953 
954  crho = self%crho ; ct = self%cT ; cb = self%cB
955 
956  eta_ohm = 0.d0 ; eta_ad = 0.d0
957  do i=1,n(1)
958  do j=1,n(2)
959  do k=1,n(3)
960  call comp_weight_element(rho(i,j,k),self%rhoarr ,crho ,rhoidx,w1)
961  call comp_weight_element(t(i,j,k) ,self%Tarr ,ct ,tidx ,w2)
962  call comp_weight_element(b(i,j,k) ,self%Barr ,cb ,bidx ,w3)
963  do iw=1,2
964  do jw=1,2
965  do kw=1,2
966  if (self%mhd_Ohm) then
967  eta_ohm(i,j,k) = eta_ohm(i,j,k) + w1(iw)*w2(jw)*w3(kw) * self%eta_Ohm_tbl(rhoidx+iw, tidx+jw, bidx+kw)
968  end if
969  if (self%mhd_AD) then
970  eta_ad(i,j,k) = eta_ad(i,j,k) + w1(iw)*w2(jw)*w3(kw) * self%eta_AD_tbl (rhoidx+iw, tidx+jw, bidx+kw)
971  end if
972  end do
973  end do
974  end do
975  end do
976  end do
977  end do
978 end subroutine eta_table
979 
980 subroutine comp_weight_element(x,xarr,cx,idx,w)
981  implicit none
982  real(8) :: x
983  real, dimension(:) :: w
984  real(8) :: x0,cx
985  real, dimension(:) :: xarr
986  real :: arr, xd, xu, idx_r
987  real(8) :: mult
988  integer(8) :: idx
989  integer :: nxarr
990 
991  nxarr = size(xarr)
992  idx = 0
993  mult = 0.d0
994  xd = 0.d0
995  xu = 0.d0
996  x0 = xarr(1)
997  mult = x / x0
998  ! log rule: log_b(x) = log_d(x)/log_d(b)
999  idx = int( log10(mult) / log10(cx) )
1000  idx_r = real(idx)
1001 
1002  xd = x0 * cx**idx_r
1003  xu = x0 * cx**(idx_r+1.d0)
1004  ! weight factors
1005  w(1) = ( x - xd ) / ( xu - xd )
1006  w(2) = 1. - w(1)
1007 
1008 end subroutine comp_weight_element
1009 
1010 subroutine comp_weight(x,xarr,cx,idx,w)
1011  implicit none
1012  real(8), dimension(:,:,:) :: x
1013  real, dimension(:,:,:,:) :: w
1014  real(8) :: x0,cx
1015  real, dimension(:) :: xarr
1016  real, allocatable, dimension(:,:,:) :: arr, xd, xu, idx_r
1017  real(8), allocatable, dimension(:,:,:) :: mult
1018  integer(8), dimension(:,:,:) :: idx
1019  integer :: n(3)
1020  integer :: nxarr
1021 
1022  n = shape(x)
1023  nxarr = size(xarr)
1024  allocate ( mult(n(1), n(2), n(3) ) )
1025  allocate ( arr(n(1), n(2), n(3) ) )
1026  allocate ( xd(n(1), n(2), n(3) ) )
1027  allocate ( xu(n(1), n(2), n(3) ) )
1028  allocate ( idx_r(n(1), n(2), n(3) ) )
1029  arr = 1.d0
1030  idx = 0
1031  mult = 0.d0
1032  xd = 0.d0
1033  xu = 0.d0
1034  x0 = xarr(1)
1035  mult = x / (x0 * arr)
1036  ! log rule: log_b(x) = log_d(x)/log_d(b)
1037  idx = int( log10(mult) / log10(cx*arr) )
1038  idx_r = real(idx)
1039 
1040  xd = x0*arr * (cx*arr)**idx_r
1041  xu = x0*arr * (cx*arr)**(idx_r+1.d0)
1042  ! weight factors
1043  w(:,:,:,1) = ( x - xd ) / ( xu - xd )
1044  w(:,:,:,2) = arr - w(:,:,:,1)
1045 
1046 
1047  deallocate(mult,arr,xd,xu,idx_r)
1048 
1049 end subroutine comp_weight
1050 
1051 subroutine eta_table_init(self,eta_ohm_tbl,eta_ad_tbl)
1052  class(non_ideal_t):: self
1053  real, dimension(:,:,:) :: eta_ohm_tbl, eta_ad_tbl
1054  real :: rho0,rhoend,crho,t0,tend,ct,e0,eend,ce,b0,bend,cb
1055  real, dimension(200) :: idxarr,rhoarr
1056  real, dimension(60) :: tarr
1057  real, dimension(150) :: barr
1058  real, dimension(60) :: earr
1059  integer :: i,j
1060  real :: nrho, nt, nb, ne, kb, mu, mp, na
1061 
1062  kb = self%kB ! Boltzmann constant
1063  mu = self%mu ! mean molecular weight
1064  mp = self%mp ! Proton mass
1065  na = self%NA ! Avogadro constant
1066 
1067  nrho = 200. ; nt = 60. ; ne = nt ; nb = 150.
1068  rho0 = 300. * mu * mp ; rhoend = 1.7688710202180409e25 * mu * mp
1069  crho= (rhoend/rho0)**(1. / (nrho - 1.) )
1070 
1071  t0 = 10. ; tend = 1.2449399509181750e5
1072  ct = (tend/t0)**(1. / (nt - 1.) )
1073 
1074  e0 = 3./2 * kb * na * t0 ; eend = 3./2 * kb * na * tend
1075  ce = (eend/e0)**(1. / (ne - 1.) )
1076 
1077  b0 = 1e-10 ; bend = 1e10
1078  cb = (bend/b0)**(1. / (nb - 1.) )
1079 
1080  do i=1,max(int(nrho),int(nt),int(nb))
1081  idxarr(i) = real(i-1)
1082  end do
1083  self%rhoarr = rho0 * crho**idxarr(1:int(nrho))
1084  self%Tarr = t0 * ct**idxarr(1:int(nt))
1085  self%Earr = e0 * ce**idxarr(1:int(ne))
1086  self%Barr = b0 * cb**idxarr(1:int(nb))
1087 
1088  self%rho0 = rho0
1089  self%crho = crho
1090  self%T0 = t0
1091  self%cT = ct
1092  self%E0 = e0
1093  self%cE = ce
1094  self%B0 = b0
1095  self%cB = cb
1096 
1097  print*,'call read table'
1098  call read_table(eta_ohm_tbl,eta_ad_tbl)
1099  print*,'read table was called!'
1100 
1101 end subroutine eta_table_init
1102 
1103 
1104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1105 ! Subroutine to read in resistivites from Marchand et al. 2016
1106 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1107 subroutine read_table(eta_ohm_tbl, eta_ad_tbl)
1108 implicit none
1109 integer :: alpha ! Number of bins
1110 integer :: nvarchimie ! Number of species
1111 integer :: nvar ! Equal to Nvarchimie
1112 integer :: nion=9 ! Number of non-grain species
1113 
1114 integer,parameter::dp=kind(1.0d0) ! default
1115 !!!! Grains info
1116 real(dp), allocatable, dimension(:) :: x_g ! Abundance of grains per bins
1117 real(dp), allocatable, dimension(:) :: r_g ! Radius of grains per bins
1118 real(dp), allocatable, dimension(:) :: m_g ! Mass of grains per bins
1119 
1120 !!!! Constants
1121 real(dp), parameter :: pi=3.1415927410125732422_dp
1122 real(dp), parameter :: mp=1.6726d-24 ! Proton mass in g
1123 real(dp), parameter :: me=9.1094d-28 ! Electron mass in g
1124 real(dp), parameter :: e=4.803204d-10 ! Electron charge in cgs
1125 real(dp), parameter :: kb = 1.3807d-16 ! Boltzmann constant erg/K
1126 real(dp), allocatable, dimension(:) :: q ! Electric charge
1127 real(dp), allocatable, dimension(:) :: m ! Mass
1128 
1129 !!!! bins of grains
1130 real(dp), parameter :: rho_s=2.3_dp ! Bulk density g.cc
1131 real(dp), parameter :: a_0=0.0375d-4 ! Reference radius cm
1132 real(dp), parameter :: a_min=0.0181d-4 ! Minimum radius cm
1133 real(dp), parameter :: a_max=0.9049d-4 ! Maximum radius cm
1134 real(dp), parameter :: zeta=a_min/a_max ! a_min/a_max
1135 real(dp), parameter :: lambda_pow=-3.5d0 ! Coeff power law
1136 real(dp) :: rho_gtot ! Grain density
1137 real(dp) :: lp1,lp3,lp4
1138 
1139 ! resistivites (cf Kunz & Mouschovias 2009, Marchand et al. 2016)
1140 real(dp), allocatable, dimension(:) :: sigma
1141 real(dp), allocatable, dimension(:) :: zetas
1142 real(dp), allocatable, dimension(:) :: phi
1143 real(dp), allocatable, dimension(:) :: tau_sn,tau_gp,tau_gm
1144 real(dp), allocatable, dimension(:,:) :: tau_inel
1145 real(dp), allocatable, dimension(:) :: gamma_zeta,gamma_omega
1146 real(dp), allocatable, dimension(:) :: omega,omega_bar
1147 real(dp), parameter :: c_l=299792458.d2 ! c in cm/s
1148 real(dp) :: sigv, muuu
1149 integer :: ib,ih,it,ix
1150 real(dp) :: b,bmaxchimie,nh,t,xi,sigh,sigo,sigp
1151 
1152 integer :: cmp_x
1153 
1154 ! For loops
1155 integer :: tchimie ! tchmie steps in temperature
1156 real(dp) :: dtchimie ! dtchmie step in temperature
1157 real(dp) :: tminchimie ! min tchmie
1158 integer :: nchimie ! nchmie steps in density
1159 real(dp) :: dnchimie ! dnchmie step in density
1160 real(dp) :: nminchimie ! min nchmie
1161 integer :: bchimie ! bchmie steps in B field
1162 real(dp) :: dbchimie ! dbchmie step in B field
1163 real(dp) :: bminchimie ! min bchmie
1164 integer :: xichimie ! Steps in ionisation rate
1165 real(dp) :: dxichimie ! dxichimie step in ion rate
1166 real(dp) :: ximinchimie ! min xichimie
1167 integer :: nislin,tislin,xiislin ! Linear scale or not
1168 real(dp),allocatable,dimension(:,:,:,:,:)::conductivities ! resistivities
1169 real(dp),allocatable,dimension(:,:,:,:)::abundances ! abundances
1170 integer :: i,j,k,ij
1171 
1172 
1173 real(dp) :: sig_p, sig_perp, sig_h, eta_ohm, eta_hall, eta_ad
1174 
1175 real, dimension(:,:,:) :: eta_ohm_tbl, eta_ad_tbl
1176 
1177  open(42,file='../../tables/non-ideal_tbl/Table_abundances.dat', status='old')
1178  read(42,*) nchimie, tchimie, xichimie, nvar ! Number of steps in density, temperature and ionisation rate
1179  read(42,*) nislin, tislin, xiislin ! 1: Constant step in logscale, 0: Not constant step in logscale
1180  read(42,*)
1181  read(42,*)
1182  read(42,*)
1183  allocate(abundances(-2:nvar,nchimie,tchimie,xichimie))
1184  ! abundances(-2,:,:,:)=density
1185  ! abundances(-1,:,:,:)=temperature
1186  ! abundances(0,:,:,:)=ionisation rate
1187  ! abundances(n>°,:,:,:)=relative abundance of species n
1188 
1189 
1190  do ij=1,xichimie
1191  do i=1,tchimie
1192  do j=1,nchimie
1193  read(42,*)abundances(-2:nvar,j,i,ij)
1194  end do
1195  read(42,*)
1196  read(42,*)
1197  end do
1198  end do
1199  close(42)
1200  nminchimie=(abundances(-2,1,1,1)) !min density
1201  tminchimie=(abundances(-1,1,1,1)) !min temperature
1202  ximinchimie=(abundances(0,1,1,1)) !min ionisation rate
1203  ! Density, temperature and IR steps in log (if nislin, tislin, xiislin = 1)
1204  dnchimie=(log10(abundances(-2,nchimie,1,1))-log10(abundances(-2,1,1,1)))/(nchimie-1)
1205  dtchimie=(log10(abundances(-1,1,tchimie,1))-log10(abundances(-1,1,1,1)))/(tchimie-1)
1206  dxichimie=(log10(abundances(0,1,1,xichimie))-log10(abundances(0,1,1,1)))/(xichimie-1)
1207 
1208 
1209  !Allocating arrays
1210  alpha=(nvar-nion)/3
1211 
1212  allocate(x_g(alpha))
1213  allocate(r_g(alpha))
1214  allocate(m_g(alpha))
1215  allocate(q(nvar))
1216  allocate(m(nvar))
1217  allocate(sigma(nvar))
1218  allocate(zetas(nvar))
1219  allocate(phi(nvar))
1220  allocate(tau_sn(nvar))
1221  allocate(tau_gp(nvar))
1222  allocate(tau_gm(nvar))
1223  allocate(gamma_zeta(nvar))
1224  allocate(gamma_omega(nvar))
1225  allocate(omega(nvar))
1226  allocate(omega_bar(nvar))
1227  allocate(tau_inel(nvar,nvar))
1228 
1229 
1230  ! Computation of grains radii and species mass
1231  lp1=dble(lambda_pow+1)
1232  lp3=dble(lambda_pow+3)
1233  lp4=dble(lambda_pow+4)
1234 
1235  ! Grain radius
1236  if (alpha==1) then
1237  r_g(1)=a_0
1238  else
1239  do i=1,alpha ! cf Marchand et al. 2016
1240  r_g(i)=a_min*zeta**(-dble(i)/dble(alpha)) * &
1241  & dsqrt( lp1/lp3* (1d0-zeta**(lp3/dble(alpha)))/(1d0-zeta**(lp1/dble(alpha))))
1242  end do
1243  end if
1244 
1245 
1246  q(:)=1.d0*e ! cations
1247  q(1)=-1.d0*e ! electrons
1248  do i=nion+1,nvar
1249  if (mod(i-nion,3)==1) q(i)=1.d0*e ! g+
1250  if (mod(i-nion,3)==2) q(i)=-1.d0*e ! g-
1251  if (mod(i-nion,3)==0) q(i)=0.d0 ! g0
1252  end do
1253  m(:)=0.d0
1254  m(1) = me ! e-
1255  m(2) = 23.5d0*mp ! ions metalliques
1256  m(3) = 29.d0*mp ! ions moleculaires
1257  m(4) = 3*mp ! H3+
1258  m(5) = mp ! H+
1259  m(6) = 12.d0*mp ! C+
1260  m(7) = 4.d0*mp ! He+
1261  m(8) = 39.098*mp ! K+
1262  m(9) = 22.99d0*mp ! Na+
1263  do i=1,alpha ! masse des grains
1264  m_g(i)=4.d0/3.d0*pi*r_g(i)**3*rho_s
1265  m(nion+3*(i-1)+1:nion+3*i)=m_g(i)
1266  end do
1267 
1268 
1269  ! values for Btable
1270  bminchimie=1.d-10
1271  bmaxchimie=1.d10 ! ok for first core in nimhd. maybe not enough for second core.
1272  bchimie=150
1273  dbchimie=(log10(bmaxchimie)-log10(bminchimie))/real((bchimie-1),dp)
1274 
1275  allocate(conductivities(0:3,1:nchimie,1:tchimie,1:xichimie,1:bchimie))
1276  !conductivities(1,:,:,:)= sigma_//
1277  !conductivities(2,:,:,:)= sigma_perp
1278  !conductivities(3,:,:,:)= sigma_hall
1279  !conductivities(0,:,:,:)= sigma_hall sign
1280 
1281 
1282  ! Computation of resistivities
1283 
1284  tau_sn = 0.0_dp
1285  omega = 0.0_dp
1286  sigma = 0.0_dp
1287  phi = 0.0_dp
1288  zetas = 0.0_dp
1289  gamma_zeta = 0.0_dp
1290  gamma_omega = 0.0_dp
1291  omega_bar = 0.0_dp
1292 
1293 
1294 do ib=1,bchimie
1295 do ix=7,7 ! 7 : cosmic-ray ionisation rate of 1e-17 s^-1
1296 !do iX=1,xichimie
1297 do it=1,tchimie
1298 do ih=1,nchimie
1299 
1300  nh=abundances(-2,ih,it,ix) ! density (.cc) of current point
1301  b =10.0d0**(log10(bminchimie)+dble(ib-1)*dbchimie) ! Magnetic field
1302  t =abundances(-1,ih,it,ix) ! Temperature
1303  xi =abundances(0,ih,it,ix) ! Ionisation rate
1304 
1305 
1306  do i=1,nion
1307  if (i==1) then ! electron
1308  sigv=3.16d-11 * (dsqrt(8d0*kb*1d-7*t/(pi*me*1d-3))*1d-3)**1.3d0
1309  tau_sn(i) = 1.d0/1.16d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nh/2.d0*sigv)
1310  else if (i>=2 .and. i<=nion) then ! ions
1311  muuu=m(i)*2d0*mp/(m(i)+2d0*mp)
1312  if (i==2 .or. i==3) then
1313  sigv=2.4d-9 *(dsqrt(8d0*kb*1d-7*t/(pi*muuu*1d-3))*1d-3)**0.6d0
1314  else if (i==4) then
1315  sigv=2d-9 * (dsqrt(8d0*kb*1d-7*t/(pi*muuu*1d-3))*1d-3)**0.15d0
1316  else if (i==5) then
1317  sigv=3.89d-9 * (dsqrt(8d0*kb*1d-7*t/(pi*muuu*1d-3))*1d-3)**(-0.02d0)
1318  else
1319  sigv=1.69d-9
1320  end if
1321  tau_sn(i) = 1.d0/1.14d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nh/2.d0*sigv)
1322  end if
1323  omega(i) = q(i)*b/(m(i)*c_l)
1324  sigma(i) = abundances(i,ih,it,ix)*nh*(q(i))**2*tau_sn(i)/m(i)
1325  phi(i) = 0.d0
1326  zetas(i) = 0.d0
1327  gamma_zeta(i) = 1.d0
1328  gamma_omega(i) = 1.d0
1329  omega_bar(i) = 0.d0
1330  end do
1331 
1332  do i=1,alpha ! grains
1333  tau_sn(nion+1+3*(i-1))=1.d0/1.28d0*(m_g(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nh/2.d0*(pi*r_g(i)**2*(8.d0*kb*t/(pi*2.d0*mp))**0.5))
1334  omega(nion+1+3*(i-1)) = q(nion+1+3*(i-1))*b/(m_g(i)*c_l)
1335  sigma(nion+1+3*(i-1)) = abundances(nion+1+3*(i-1),ih,it,ix)*nh*(q(nion+1+3*(i-1)))**2*tau_sn(nion+1+3*(i-1))/m_g(i)
1336 
1337  tau_sn(nion+2+3*(i-1))=tau_sn(nion+1+3*(i-1))
1338  omega(nion+2+3*(i-1)) = q(nion+2+3*(i-1))*b/(m_g(i)*c_l)
1339  sigma(nion+2+3*(i-1)) = abundances(nion+2+3*(i-1),ih,it,ix)*nh*(q(nion+2+3*(i-1)))**2*tau_sn(nion+2+3*(i-1))/m_g(i)
1340 
1341  end do
1342 
1343  sigp=0.d0
1344  sigo=0.d0
1345  sigh=0.d0
1346 
1347  do i=1,nvar
1348  sigp=sigp+sigma(i)
1349  sigo=sigo+sigma(i)/(1.d0+(omega(i)*tau_sn(i))**2)
1350  sigh=sigh-sigma(i)*omega(i)*tau_sn(i)/(1.d0+(omega(i)*tau_sn(i))**2)
1351  end do
1352 
1353  conductivities(1,ih,it,ix,ib)=log10(sigp)
1354  conductivities(2,ih,it,ix,ib)=log10(sigo)
1355  !conductivities(3,iH,iT,iX,iB)=log10(abs(sigH))
1356  !conductivities(0,iH,iT,iX,iB)=sign(1.0_dp,sigH)
1357 
1358  !conductivities(:,:,:,:) contains the log10(conductivities)
1359 
1360 
1361  eta_ohm = 1d0/sigp
1362  eta_hall = sigh/(sigo**2 + sigh**2)
1363  eta_ad = sigo/(sigo**2 + sigh**2)-1d0/sigp
1364 
1365  eta_ohm_tbl(ih,it,ib) = eta_ohm
1366  eta_ad_tbl(ih,it,ib) = eta_ad
1367  if (eta_ohm .lt. 0) eta_ohm_tbl(ih,it,ib) = 0d0
1368  if (eta_ad .lt. 0) eta_ad_tbl(ih,it,ib) = 0d0
1369 
1370 end do
1371 end do
1372 end do
1373 end do
1374 !open(unit=3,file='eta_AD_tbl.dat',ACCESS='APPEND')
1375 !write(3,*) eta_ad_tbl
1376 !close(3)
1377 
1378 end subroutine
1379 
1380 end module
1381 ! modif nimhd
1382 
1383 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
non-ideal MHD module
Template module for mesh.
Definition: mesh_mod.f90:4
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4