18 real:: eta_ohm, etamd, gamma_ad, rho_ions
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
41 procedure:: non_ideal_comp
43 procedure:: non_ideal_emf_up
44 procedure:: computejb2
45 procedure:: computambip
46 procedure:: computdifmag
48 procedure:: eta_table_init
58 SUBROUTINE init (self)
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
70 read (io_unit%input, non_ideal_params,iostat=iostat)
71 if (io%master)
write (*, non_ideal_params)
74 self%mhd_Ohm = mhd_ohm
75 self%eta_Ohm = eta_ohm
77 self%gamma_AD = gamma_ad
78 self%rho_ions = rho_ions
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
87 self%mp = 1.6737236e-24
88 self%NA = 6.022140857e23
89 call eta_table_init(self,self%eta_Ohm_tbl,self%eta_AD_tbl)
94 subroutine non_ideal_comp(self,mesh,idx,bf,uin,qin,fluxambdiff,emfambdiff,u_max,ds,dt,gamma)
98 integer:: n(4), l(3), u(3)
102 real,
dimension(:,:,:,:):: bf, uin, qin, fluxambdiff, emfambdiff
104 real(8),
allocatable,
dimension(:,:,:)::emfx
105 real(8),
allocatable,
dimension(:,:,:)::emfy
106 real(8),
allocatable,
dimension(:,:,:)::emfz
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
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))
161 if((self%mhd_AD).or.(self%mhd_Ohm))
then 163 call computejb2(self,bf,qin,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,ds)
170 call computambip(self,uin,qin,bemfx,bemfy,bemfz,florentzx,florentzy,florentzz,fluxad,bmagij,emfambdiff,fluxambdiff,jxbsquare,ds,dt,gamma)
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)) ) )
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))
189 deallocate(emfx,emfy,emfz,bmagij,jcentersquare,jxbsquare,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,emfohmdiss,fluxohm)
193 subroutine flux_upd(self,flux,fluxambdiff, dx, dt, idim)
197 integer :: i, j, k, idim
198 real,
dimension(:,:,:,:,:):: flux
199 real,
dimension(:,:,:,:)::fluxambdiff
201 real,
allocatable,
dimension(:,:,:,:)::fluxohm
206 allocate (fluxohm(n(1),n(2),n(3),3))
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
218 subroutine non_ideal_emf_up(self,emfup,emfambdiff, dx, dt, idim)
222 real,
dimension(:,:,:)::emfup
223 real,
dimension(:,:,:,:)::emfambdiff
224 real,
allocatable,
dimension(:,:,:,:)::emfohmdiss
227 n = shape(emfambdiff)
228 allocate (emfohmdiss(n(1),n(2),n(3),n(4)))
231 if (self%mhd_AD .eqv. .true. .or. self%mhd_Ohm .eqv. .true.)
then 235 emfup(:,:,:) = ( emfambdiff(:,:,:,idim) + emfohmdiss(:,:,:,idim) ) * dt / dx
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
242 deallocate(emfohmdiss)
247 subroutine computejb2(self,bf,q,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,florentzx,florentzy,florentzz,fluxmd,fluxh,fluxad,ds)
250 integer :: i, j, k, ibegin, iend, jbegin, jend, kbegin, kend
251 real,
dimension(:,:,:,:)::bf
252 real,
dimension(:,:,:,:)::q
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
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
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))
288 dx = ds(1) ; dy = ds(2) ; dz = ds(3)
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) )
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)
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) )
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) )
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) )
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)
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) )
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) )
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) )
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) )
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) )
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
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)
364 if((self%mhd_AD).or.(self%mhd_Ohm))
then 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
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,:)
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)
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,:)
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)
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,:)
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)
410 deallocate (bmagijbis, jface, bcenter, fluxbis, fluxter, fluxquat)
414 end subroutine computejb2
419 subroutine computdifmag(self,u,q,bemfx,bemfy,bemfz,jemfx,jemfy,jemfz,bmagij,fluxmd,emfohmdiss,fluxohm,jcentersquare,ds,dt)
425 real,
dimension(:,:,:,:)::u
426 real,
dimension(:,:,:,:)::q
428 real(8),
dimension(:,:,:,:)::bemfx,bemfy,bemfz
429 real(8),
dimension(:,:,:,:)::jemfx,jemfy,jemfz
430 real(8),
dimension(:,:,:,:,:)::bmagij
431 real(8),
dimension(:,:,:,:)::fluxmd
433 real(8),
dimension(:,:,:,:):: emfohmdiss,fluxohm
434 real(8),
dimension(:,:,:)::jcentersquare
442 real(8),
allocatable,
dimension(:,:,:,:)::jcenter
443 real(8),
allocatable,
dimension(:,:,:,:,:)::jface
444 real(8),
allocatable,
dimension(:,:,:,:)::jemf
446 real(8)::rhocell,bcell,cv
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
452 integer,
dimension(3) :: index_i,index_j,index_k
455 real(8):: dx, dy, dz, dt
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) )
498 dx = ds(1); dy = ds(2) ; dz = ds(3)
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))
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
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,:)
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
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))
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
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)
544 deallocate(jcenter, jface, jemf, bsquarex, bsquarey, bsquarez, dtlim, rhox, rhoy, rhoz, epsx, epsy, epsz, rhof, pf, bsqf, epsf )
546 end subroutine computdifmag
551 SUBROUTINE computambip(self,u,q,bemfx,bemfy,bemfz,florentzx,florentzy,florentzz,fluxad,bmagij,emfambdiff,fluxambdiff,jxbsquare,ds,dt,gamma)
554 real,
dimension(:,:,:,:)::u
555 real,
dimension(:,:,:,:)::q
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
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
577 real(8),
allocatable,
dimension(:,:,:,:)::florentz
578 real(8),
allocatable,
dimension(:,:,:) ::bsquaremax
580 real(8),
allocatable,
dimension(:,:,:,:)::jcenter
581 real(8),
allocatable,
dimension(:,:,:,:)::jxb
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
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
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))
673 dx = ds(1) ; dy = ds(2) ; dz = ds(3)
675 if ( self%ntestDADM .eqv. .false. )
then 679 scale_t = scale_l/scale_v
680 scale_t2 = (scale_v)**2 * self%mp / self%kB
682 c2_4pi = 7.152426325648638e+19
683 bcode_to_gauss = (4.d0*pi*scale_d)**0.5*scale_v
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))
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))
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)
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
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
712 bsquaremax(2:n(1)-1,2:n(2)-1,2:n(3)-1)=max(bsquarex,bsquarey,bsquarez,bsquarexx,bsquareyy,bsquarezz)
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)
719 rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
721 if (self%ntestDADM)
then 722 betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhox)
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 ) - &
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)
742 eta_ad = eta_ad * c2_4pi
743 eta_ad = eta_ad / scale_v**2 / scale_t
744 betaad2 = eta_ad / brms2
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
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)
754 rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
757 if (self%ntestDADM)
then 758 betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhoy)
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 ) - &
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)
777 eta_ad = eta_ad * c2_4pi
778 eta_ad = eta_ad / scale_v**2 / scale_t
779 betaad2 = eta_ad / brms2
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
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)
788 rhocell = rhocellmin(2:n(1)-1,2:n(2)-1,2:n(3)-1)
790 if (self%ntestDADM)
then 791 betaad2=1.d0/(self%gamma_AD*self%rho_ions*rhoz)
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 ) - &
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)
810 eta_ad = eta_ad * c2_4pi
811 eta_ad = eta_ad / scale_v**2 / scale_t
812 betaad2 = eta_ad / brms2
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
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)
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))
829 eth = etfx - 0.5 / rhofx * ( uxfx**2 + uyfx**2 + uzfx**2 ) - &
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)
839 eta_ad = eta_ad * c2_4pi
840 eta_ad = eta_ad / scale_v**2 / scale_t
841 betaad2 = eta_ad / brms2
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)
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)
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))
854 eth = etfy - 0.5 / rhofy * ( uxfy**2 + uyfy**2 + uzfy**2 ) - &
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)
864 eta_ad = eta_ad * c2_4pi
865 eta_ad = eta_ad / scale_v**2 / scale_t
866 betaad2 = eta_ad / brms2
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)
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)
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))
879 eth = etfz - 0.5 / rhofz * ( uxfz**2 + uyfz**2 + uzfz**2 ) - &
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)
889 eta_ad = eta_ad * c2_4pi
890 eta_ad = eta_ad / scale_v**2 / scale_t
891 betaad2 = eta_ad / brms2
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)
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
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)
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))
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 ) - &
913 rhocell = u(2:n(1)-1,2:n(2)-1,2:n(3)-1,1)
916 rho = rhocell * scale_d
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)
923 eta_ad = eta_ad * c2_4pi
924 eta_ad = eta_ad / scale_v**2 / scale_t
925 betaad = eta_ad / brms2
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))*&
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 )
938 end SUBROUTINE computambip
942 subroutine eta_table(self,eta_AD,eta_Ohm,rho,T,B)
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
954 crho = self%crho ; ct = self%cT ; cb = self%cB
956 eta_ohm = 0.d0 ; eta_ad = 0.d0
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)
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)
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)
978 end subroutine eta_table
980 subroutine comp_weight_element(x,xarr,cx,idx,w)
983 real,
dimension(:) :: w
985 real,
dimension(:) :: xarr
986 real :: arr, xd, xu, idx_r
999 idx = int( log10(mult) / log10(cx) )
1003 xu = x0 * cx**(idx_r+1.d0)
1005 w(1) = ( x - xd ) / ( xu - xd )
1008 end subroutine comp_weight_element
1010 subroutine comp_weight(x,xarr,cx,idx,w)
1012 real(8),
dimension(:,:,:) :: x
1013 real,
dimension(:,:,:,:) :: w
1015 real,
dimension(:) :: xarr
1016 real,
allocatable,
dimension(:,:,:) :: arr, xd, xu, idx_r
1017 real(8),
allocatable,
dimension(:,:,:) :: mult
1018 integer(8),
dimension(:,:,:) :: idx
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) ) )
1035 mult = x / (x0 * arr)
1037 idx = int( log10(mult) / log10(cx*arr) )
1040 xd = x0*arr * (cx*arr)**idx_r
1041 xu = x0*arr * (cx*arr)**(idx_r+1.d0)
1043 w(:,:,:,1) = ( x - xd ) / ( xu - xd )
1044 w(:,:,:,2) = arr - w(:,:,:,1)
1047 deallocate(mult,arr,xd,xu,idx_r)
1049 end subroutine comp_weight
1051 subroutine eta_table_init(self,eta_ohm_tbl,eta_ad_tbl)
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
1060 real :: nrho, nt, nb, ne, kb, mu, mp, na
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.) )
1071 t0 = 10. ; tend = 1.2449399509181750e5
1072 ct = (tend/t0)**(1. / (nt - 1.) )
1074 e0 = 3./2 * kb * na * t0 ; eend = 3./2 * kb * na * tend
1075 ce = (eend/e0)**(1. / (ne - 1.) )
1077 b0 = 1e-10 ; bend = 1e10
1078 cb = (bend/b0)**(1. / (nb - 1.) )
1080 do i=1,max(int(nrho),int(nt),int(nb))
1081 idxarr(i) =
real(i-1)
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))
1097 print*,
'call read table' 1098 call read_table(eta_ohm_tbl,eta_ad_tbl)
1099 print*,
'read table was called!' 1101 end subroutine eta_table_init
1107 subroutine read_table(eta_ohm_tbl, eta_ad_tbl)
1110 integer :: nvarchimie
1114 integer,
parameter::dp=kind(1.0d0)
1116 real(dp),
allocatable,
dimension(:) :: x_g
1117 real(dp),
allocatable,
dimension(:) :: r_g
1118 real(dp),
allocatable,
dimension(:) :: m_g
1121 real(dp),
parameter :: pi=3.1415927410125732422_dp
1122 real(dp),
parameter :: mp=1.6726d-24
1123 real(dp),
parameter :: me=9.1094d-28
1124 real(dp),
parameter :: e=4.803204d-10
1125 real(dp),
parameter :: kb = 1.3807d-16
1126 real(dp),
allocatable,
dimension(:) :: q
1127 real(dp),
allocatable,
dimension(:) :: m
1130 real(dp),
parameter :: rho_s=2.3_dp
1131 real(dp),
parameter :: a_0=0.0375d-4
1132 real(dp),
parameter :: a_min=0.0181d-4
1133 real(dp),
parameter :: a_max=0.9049d-4
1134 real(dp),
parameter :: zeta=a_min/a_max
1135 real(dp),
parameter :: lambda_pow=-3.5d0
1136 real(dp) :: rho_gtot
1137 real(dp) :: lp1,lp3,lp4
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
1148 real(dp) :: sigv, muuu
1149 integer :: ib,ih,it,ix
1150 real(dp) :: b,bmaxchimie,nh,t,xi,sigh,sigo,sigp
1156 real(dp) :: dtchimie
1157 real(dp) :: tminchimie
1159 real(dp) :: dnchimie
1160 real(dp) :: nminchimie
1162 real(dp) :: dbchimie
1163 real(dp) :: bminchimie
1165 real(dp) :: dxichimie
1166 real(dp) :: ximinchimie
1167 integer :: nislin,tislin,xiislin
1168 real(dp),
allocatable,
dimension(:,:,:,:,:)::conductivities
1169 real(dp),
allocatable,
dimension(:,:,:,:)::abundances
1173 real(dp) :: sig_p, sig_perp, sig_h, eta_ohm, eta_hall, eta_ad
1175 real,
dimension(:,:,:) :: eta_ohm_tbl, eta_ad_tbl
1177 open(42,file=
'../../tables/non-ideal_tbl/Table_abundances.dat', status=
'old')
1178 read(42,*) nchimie, tchimie, xichimie, nvar
1179 read(42,*) nislin, tislin, xiislin
1183 allocate(abundances(-2:nvar,nchimie,tchimie,xichimie))
1193 read(42,*)abundances(-2:nvar,j,i,ij)
1200 nminchimie=(abundances(-2,1,1,1))
1201 tminchimie=(abundances(-1,1,1,1))
1202 ximinchimie=(abundances(0,1,1,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)
1212 allocate(x_g(alpha))
1213 allocate(r_g(alpha))
1214 allocate(m_g(alpha))
1217 allocate(sigma(nvar))
1218 allocate(zetas(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))
1231 lp1=dble(lambda_pow+1)
1232 lp3=dble(lambda_pow+3)
1233 lp4=dble(lambda_pow+4)
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))))
1249 if (mod(i-nion,3)==1) q(i)=1.d0*e
1250 if (mod(i-nion,3)==2) q(i)=-1.d0*e
1251 if (mod(i-nion,3)==0) q(i)=0.d0
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)
1273 dbchimie=(log10(bmaxchimie)-log10(bminchimie))/
real((bchimie-1),dp)
1275 allocate(conductivities(0:3,1:nchimie,1:tchimie,1:xichimie,1:bchimie))
1290 gamma_omega = 0.0_dp
1300 nh=abundances(-2,ih,it,ix)
1301 b =10.0d0**(log10(bminchimie)+dble(ib-1)*dbchimie)
1302 t =abundances(-1,ih,it,ix)
1303 xi =abundances(0,ih,it,ix)
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 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
1315 sigv=2d-9 * (dsqrt(8d0*kb*1d-7*t/(pi*muuu*1d-3))*1d-3)**0.15d0
1317 sigv=3.89d-9 * (dsqrt(8d0*kb*1d-7*t/(pi*muuu*1d-3))*1d-3)**(-0.02d0)
1321 tau_sn(i) = 1.d0/1.14d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nh/2.d0*sigv)
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)
1327 gamma_zeta(i) = 1.d0
1328 gamma_omega(i) = 1.d0
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)
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)
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)
1353 conductivities(1,ih,it,ix,ib)=log10(sigp)
1354 conductivities(2,ih,it,ix,ib)=log10(sigo)
1362 eta_hall = sigh/(sigo**2 + sigh**2)
1363 eta_ad = sigo/(sigo**2 + sigh**2)-1d0/sigp
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
Template module for mesh.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Module with list handling for generic class task_t objects.