8 subroutine calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
24 integer,
intent(in) :: qunit, igrid
25 double precision,
intent(in),
dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
26 double precision,
intent(in),
dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
27 integer :: ixC^L,ixCC^L
28 logical,
intent(in) :: first
30 double precision,
dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
31 double precision,
dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
32 double precision,
dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
33 double precision,
dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
34 double precision,
dimension(0:nw+nwauxio),
intent(out) :: normconv
36 double precision :: ldw(ixG^T), dwC(ixG^T)
37 double precision,
dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
38 double precision,
dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
39 double precision,
dimension(ixG^T,1:nw+nwauxio) :: w
40 integer :: nxCC^D,idims,jxC^L,iwe
41 integer :: nx^D, nxC^D, ix^D, ix, iw, level, idir
42 logical,
save :: subfirst=.true.
47 ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
48 ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
50 nx^d=ixmhi^d-ixmlo^d+1;
56 w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
66 do ix^d=jxcmin^d,jxcmax^d
67 w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmin^d-1^d%jxC^s,nw-nwextra+1:nw)
71 do ix^d=jxcmin^d,jxcmax^d
72 w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmax^d+1^d%jxC^s,nw-nwextra+1:nw)
88 call mpistop(
"usr_aux_output not defined")
98 if(
b0field .and.
allocated(iw_mag))
then
101 w(ixg^t,iw_e)=w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=
ndim+1) &
102 + sum(w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=
ndim+1)
104 w(ixg^t,iw_mag(:))=w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
108 wcc(ixcc^s,:)=w(ixcc^s,:)
115 {
do ix^db=ixcmin^db,ixcmax^db\}
116 wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw))/dble(2**
ndim)
121 {
do ix^db=ixcmin^db,ixcmax^db\}
122 wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw)*ps(igrid)%dvolume(ix^d:ix^d+1)) &
123 /sum(ps(igrid)%dvolume(ix^d:ix^d+1))
130 xc_tmp(ixc^s,1:ndim) = xc(ixc^s,1:ndim)
131 wc_tmp(ixc^s,1:nw+nwauxio) = wc(ixc^s,1:nw+nwauxio)
132 xcc_tmp(ixcc^s,1:ndim) = xcc(ixcc^s,1:ndim)
133 wcc_tmp(ixcc^s,1:nw+nwauxio) = wcc(ixcc^s,1:nw+nwauxio)
143 if(nwaux>0 .and. mype==0 .and. first.and.subfirst)
then
145 if(convert_type==
'idl'.or.convert_type==
'tecplot' &
146 .or.convert_type==
'vtu'.or.convert_type==
'vtuB') &
147 write(*,*)
'Warning: also averaged auxiliaries within calc_grid'
161 double precision,
dimension(ix^S,ndim) :: xC
162 double precision,
dimension(ix^S,nw+nwauxio) :: wC
163 double precision,
dimension(ix^S,ndim) :: x_TMP
164 double precision,
dimension(ix^S,nw+nwauxio) :: w_TMP
167 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
168 double precision,
dimension(ndim,ndim) :: normal
169 integer,
dimension(nw) :: vectoriw
170 integer :: ix^D, idim, iw, ivector, iw0
177 vectoriw(iw_vector(ivector)+idim)=iw_vector(ivector)
182 {
do ix^db=ixmin^db,ixmax^db\}
185 x_tec(1:ndim)=xc(ix^d,1:ndim)
186 w_tec(1:nw+nwauxio)=wc(ix^d,1:nw+nwauxio)
193 x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,2))
194 x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))
200 x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,phi_))
201 x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,phi_))
202 x_tec(3)=xc(ix^d,z_)}
205 {^ifoned normal(1,1)=one}
210 normal(1,1)=cos(xc(ix^d,2))
211 normal(1,2)=-sin(xc(ix^d,2))
212 normal(2,1)=sin(xc(ix^d,2))
213 normal(2,2)=cos(xc(ix^d,2))
222 normal(1,1)=cos(xc(ix^d,phi_))
223 normal(1,phi_)=-sin(xc(ix^d,phi_))
226 normal(2,1)=sin(xc(ix^d,phi_))
227 normal(2,phi_)=cos(xc(ix^d,phi_))
235 if (iw<=nw) iw0=vectoriw(iw)
236 if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw)
then
238 w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
240 w_tec(iw)=wc(ix^d,iw)
244 x_tec(1)=xc(ix^d,1){^nooned*sin(xc(ix^d,2))}{^ifthreed*cos(xc(ix^d,3))}
246 x_tec(2)=xc(ix^d,1)*cos(xc(ix^d,2))}
248 x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))*sin(xc(ix^d,3))
249 x_tec(3)=xc(ix^d,1)*cos(xc(ix^d,2))}
252 {^ifoned normal(1,1)=one}
254 normal(1,1)=sin(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))}
255 normal(1,2)=cos(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))
256 normal(1,3)=-sin(xc(ix^d,3))}}
259 normal(2,1)=cos(xc(ix^d,2))
260 normal(2,2)=-sin(xc(ix^d,2))}
262 normal(2,1)=sin(xc(ix^d,2))*sin(xc(ix^d,3))
263 normal(2,2)=cos(xc(ix^d,2))*sin(xc(ix^d,3))
264 normal(2,3)=cos(xc(ix^d,3))
266 normal(3,1)=cos(xc(ix^d,2))
267 normal(3,2)=-sin(xc(ix^d,2))
271 if (iw<=nw) iw0=vectoriw(iw)
272 if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw)
then
274 w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
276 w_tec(iw)=wc(ix^d,iw)
280 write(*,*)
"No converter for coordinate=",coordinate
282 x_tmp(ix^d,1:ndim)=x_tec(1:ndim)
283 w_tmp(ix^d,1:nw+nwauxio)=w_tec(1:nw+nwauxio)
285 where(dabs(w_tmp(ix^d,1:nw+nwauxio))<smalldouble)
286 w_tmp(ix^d,1:nw+nwauxio)=zero
304 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
305 character(len=1024) :: outfilehead
307 integer:: space_position,iw,ind
308 logical,
save:: first=.true.
309 character(len=name_len):: wname
310 character(len=std_len):: aux_variable_names
311 character(len=std_len):: scanstring
316 call mpistop(
"usr_add_aux_names not defined")
324 scanstring=trim(aux_variable_names)
325 wnamei(1:nw)=prim_wnames(1:nw)
327 scanstring=trim(aux_variable_names)
328 wnamei(1:nw)=cons_wnames(1:nw)
331 space_position=index(scanstring,
' ')
332 do iw=nw+1,nw+nwauxio
333 do while (space_position==1)
334 scanstring=scanstring(2:)
335 space_position=index(scanstring,
' ')
337 wname=scanstring(:space_position-1)
338 scanstring=scanstring(space_position+1:)
339 space_position=index(scanstring,
' ')
342 wnamei(iw)=trim(wname)
348 xandwnamei(1)=
"r";{^nooned xandwnamei(2)=
"Theta"};{^ifthreed xandwnamei(3)=
"Phi"}
364 xandwnamei(1)=
"X";{^nooned xandwnamei(2)=
"Y"};{^ifthreed xandwnamei(3)=
"Z"}
367 xandwnamei(ndim+1:ndim+nw+nwauxio)=wnamei(1:nw+nwauxio)
371 write(outfilehead,
'(a)') trim(xandwnamei(1))
375 write(outfilehead,
'(a)')outfilehead(1:len_trim(outfilehead))//
" "//trim(wname)
382 write(outfilehead,
'(a)')outfilehead(1:len_trim(outfilehead))//
" "//trim(wname)
387 do iw=ndim+nw+1,ndim+nw+nwauxio
389 write(outfilehead,
'(a)')outfilehead(1:len_trim(outfilehead))//
" "//trim(wname)
393 if(first.and.
mype==0)
then
394 print*,
'-------------------------------------------------------------------------------'
395 write(
unitterm,*)
'Saving visual data. Coordinate directions and variable names are:'
399 print *,ind,xandwnamei(iw)
404 write(*,*) ind,wnamei(iw-ndim)
407 do iw=ndim+nw+1,ndim+nw+nwauxio
409 print *,ind,wnamei(iw-ndim)
412 print*,
'-------------------------------------------------------------------------------'