J1xBq=fltarr(nx,ny,nz,3,m)
;
for i=0,m-1 do begin
  J1xBq(*,*,*,*,i)=cross(J1,Bi(*,*,*,*,i))
endfor
;
TM_pq=fltarr(m,m)
for i=0,m-1 do begin
for j=0,m-1 do begin
  TM_pq(i,j)=-mean(dot(ui(*,*,*,*,i),J1xBq(*,*,*,*,j)))
  print,i,j
endfor
endfor
;
TM_p=total(TM_pq,2)
END