mod_trace_field.t

Go to the documentation of this file.

module mod_trace_field
  use mod_global_parameters
  use mod_physics
  implicit none
 
contains
 
  subroutine trace_field_multi(xfm,wPm,wLm,dL,numL,numP,nwP,nwL,forwardm,ftype,tcondi)
    ! trace multiple field lines
    ! xfm: locations of points at the field lines. User should provide xfm(1:numL,1,1:ndim)
    !   as seed points, then field line wills be traced from the seed points. xfm(1:numL,1,:) 
    !   are user given and other points are given by the subroutine as feedback.
    ! numL: number of field lines user wants to trace. user given
    ! numP: maximum number of points at the field line. User defined. Note that not every
    !   point of numP is valid, as the tracing will stop when the field line leave the
    !   simulation box. The number of valid point is given in wLm(1)
    ! wPm: point variables, which have values at each point of xfm. The way to get wP is
    !   user defined, with help of IO given by subroutine set_field_w in mod_usr_methods.
    !   User can calculate the density/temperature at the point and then store the values
    !   into wPm, or do something else.
    ! nwP: number of point variables. user given
    ! wLm: line variables, variables for the lines rather than for each point of the lines.
    !   For example, wLm(1) stores the number of valid points at the field lines. The way to 
    !   get wLm is also user defined with set_field_w, the same as wPm. User can calculate
    !   the maximum density/temperature at the field lines and stores it in wLm
    ! nwL: number of line variables. user given
    ! dL: length step for field line tracing. user given
    ! forwardm: true--trace field forward; false--trace field line backward. user given
    ! ftype: type of field user wants to trace. User can trace velocity field by setting
    !   ftype='Vfield' or trace magnetic field by setting ftype='Bfield'. It is possible 
    !   to trace other fields, e.g. electric field, where user can define the field with
    !   IO given by subroutine set_field in mod_usr_methods. user given
    ! tcondi: user given
    use mod_particle_base
 
    integer, intent(in) :: numL,numP,nwP,nwL
    double precision, intent(inout) :: xfm(numL,numP,ndim),wPm(numL,numP,nwP),wLm(numL,1+nwL)
    double precision, intent(in) :: dL
    logical, intent(in) :: forwardm(numL)
    character(len=std_len), intent(in) :: ftype,tcondi
 
    integer :: indomain,ipe_now,igrid_now,igrid,j,iL
    integer :: ipoint_in,ipoint_out,numSend,nRT,nRTmax
    double precision :: x3d(3),statusF(4+ndim),statusL(numL,4+ndim+nwL),statusS(numL,4+ndim+nwL)
    logical :: continueL(numL),myL(numL)
    logical :: stopT,forward
    integer :: ipointm(numL),igridm(numL)
    double precision :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
    double precision, allocatable :: data_send(:,:,:),data_recv(:,:,:)
 
    if (tcondi/='TRAC') then
      wpm=zero
    else
      wpm=-1
    endif
    wlm=zero
    xfm(1:numl,2:nump,:)=zero
    stopt=.true.
    myl=.false.
    xf=zero
    wp=zero
 
    ! find the pe and igrid for the first point
    do il=1,numl
      indomain=0
      wlm(il,1)=0
      {if (xfm(il,1,^db)>=xprobmin^db .and. xfm(il,1,^db)<xprobmax^db) indomain=indomain+1\}
      if (indomain==ndim) then
        if (tcondi/='TRAC') wlm(il,1)=1
        continuel(il)=.true.
        ! find pe and igrid
        x3d=0.d0
        do j=1,ndim
          x3d(j)=xfm(il,1,j)
        enddo
        call find_particle_ipe(x3d,igrid_now,ipe_now)
        stopt=.false.
        ipointm(il)=1
        igridm(il)=igrid_now
        if (mype==ipe_now) then 
          myl(il)=.true.
        else
          xfm(il,1,:)=zero
        endif
      else
        continuel(il)=.false.
        wlm(il,1)=zero
      endif
    enddo
 
    do while(stopt .eqv. .false.)
      ! tracing multiple field lines inside pe
      statuss=zero
      do il=1,numl
        if (myl(il) .and. continuel(il)) then
          igrid=igridm(il)
          ipoint_in=ipointm(il)
          xf(ipoint_in,:)=xfm(il,ipoint_in,:)
          wl(:)=wlm(il,:)
          forward=forwardm(il)
          statusf=zero
          call find_points_in_pe(igrid,ipoint_in,xf,wp,wl,dl,nump,nwp,nwl,forward,ftype,tcondi,statusf)
          ipoint_out=int(statusf(1))
          xfm(il,ipoint_in:ipoint_out-1,:)=xf(ipoint_in:ipoint_out-1,:)
          wpm(il,ipoint_in:ipoint_out-1,:)=wp(ipoint_in:ipoint_out-1,:)
          ! status for each field line
          ! 1: index of next point
          ! 2: ipe of next point
          ! 3: igrid of next point
          ! 4: 1-> continue tracing; 0-> stop tracing
          ! 5:4+ndim: coordinate of next point
          ! 4+ndim+1:4+ndim+nwL: wL(2:1+nwL)
          ! for TRAC nwL=2 -> wL(2): current Tcoff; wL(3): Tmax
          ! for TRAC nwP=2 -> wP(:,1):ipe; wP(:,2):igrid
          statuss(il,1:4+ndim)=statusf(1:4+ndim)
          statuss(il,4+ndim+1:4+ndim+nwl)=wl(2:1+nwl)
          if (tcondi=='TRAC') wlm(il,1)=ipoint_out-1
        endif
      enddo
 
      ! comunicating tracing results
      numsend=numl*(4+ndim+nwl)
      call mpi_allreduce(statuss,statusl,numsend,mpi_double_precision,&
                       mpi_sum,icomm,ierrmpi)
 
      ! for next step
      stopt=.true.
      myl=.false.
      do il=1,numl
        if (continuel(il)) then
          ipointm(il)=int(statusl(il,1))
          if (mype==int(statusl(il,2))) myl(il)=.true.
          igridm(il)=int(statusl(il,3))
          if (int(statusl(il,4))==0) then
            stopt=.false.
          else
            continuel(il)=.false.
          endif
          if (myl(il)) xfm(il,ipointm(il),1:ndim)=statusl(il,4+1:4+ndim)
          if (tcondi/='TRAC') wlm(il,1)=ipointm(il)-1
          wlm(il,2:1+nwl)=statusl(il,4+ndim+1:4+ndim+nwl)
        endif
      enddo
    enddo
 
    ! communication after tracing
    if (tcondi/='TRAC') then
      nrtmax=0
      do il=1,numl
        if (nrtmax<int(wlm(il,1))) nrtmax=int(wlm(il,1))
      enddo
      numsend=numl*nrtmax*(ndim+nwp)
 
      allocate(data_send(numl,nrtmax,ndim+nwp),data_recv(numl,nrtmax,ndim+nwp))
      data_send(:,:,:)=zero
      do il=1,numl
        nrt=int(wlm(il,1))
        data_send(il,1:nrt,1:ndim)=xfm(il,1:nrt,1:ndim)
        if (nwp>0) data_send(il,1:nrt,1+ndim:ndim+nwp)=wpm(il,1:nrt,1:nwp)
      enddo
      call mpi_allreduce(data_send,data_recv,numsend,mpi_double_precision,&
                         mpi_sum,icomm,ierrmpi)
      do il=1,numl
        nrt=int(wlm(il,1))
        xfm(il,1:nrt,1:ndim)=data_recv(il,1:nrt,1:ndim)
        if (nwp>0) wpm(il,1:nrt,1:nwp)=data_recv(il,1:nrt,1+ndim:ndim+nwp)
      enddo
      deallocate(data_send,data_recv)
    endif
 
  end subroutine trace_field_multi
 
  subroutine trace_field_single(xf,wP,wL,dL,numP,nwP,nwL,forward,ftype,tcondi)
    ! trace a field line
    ! xf: locations of points at the field line. User should provide xf(1,1:ndim)
    !   as seed point, then field line will be traced from the seed point. xf(1,:) is 
    !   user given and xf(2:wL(1),:) are given by the subroutine as feedback.
    ! numP: maximum number of points at the field line. User defined. Note that not every
    !   point of numP is valid, as the tracing will stop when the field line leave the
    !   simulation box. The number of valid point is given in wL(1)
    ! wP: point variables, which have values at each point of xf. The way to get wP is
    !   user defined, with help of IO given by subroutine set_field_w in mod_usr_methods.
    !   User can calculate the density/temperature at the point and then store the values
    !   into wP, or do something else.
    ! nwP: number of point variables. user given
    ! wL: line variables, variables for the line rather than for each point of the line.
    !   For example, wL(1) stores the number of valid points at the field line. The way to 
    !   get wL is also user defined with set_field_w, the same as wP. User can calculate
    !   the maximum density/temperature at the field line and stores it in wL
    ! nwL: number of line variables. user given
    ! dL: length step for field line tracing. user given
    ! forward: true--trace field forward; false--trace field line backward. user given
    ! ftype: type of field user wants to trace. User can trace velocity field by setting
    !   ftype='Vfield' or trace magnetic field by setting ftype='Bfield'. It is possible 
    !   to trace other fields, e.g. electric field, where user can define the field with
    !   IO given by subroutine set_field in mod_usr_methods. user given
    ! tcondi: user given
    use mod_usr_methods
    use mod_particle_base
 
    integer, intent(in) :: numP,nwP,nwL
    double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
    double precision, intent(in) :: dL
    logical, intent(in) :: forward
    character(len=std_len), intent(in) :: ftype,tcondi
 
    integer :: indomain,ipoint_in,ipe_now,igrid_now,igrid,j
    integer :: ipoint_out,ipe_next,igrid_next,numRT
    double precision :: x3d(3),statusF(4+ndim),status_bcast(4+ndim+nwL)
    logical :: stopT
    double precision, allocatable :: data_send(:,:),data_recv(:,:)
 
    wp=zero
    wl=zero
    xf(2:nump,:)=zero
 
    ! check whether or the first point is inside simulation box. if yes, find
    ! the pe and igrid for the point
    indomain=0
    wl(1)=0
    {if (xf(1,^db)>=xprobmin^db .and. xf(1,^db)<xprobmax^db) indomain=indomain+1\}
    if (indomain==ndim) then
      wl(1)=1
 
       ! find pe and igrid
       x3d=0.d0
       do j=1,ndim
         x3d(j)=xf(1,j)
       enddo
      call find_particle_ipe(x3d,igrid_now,ipe_now)
      stopt=.false.
      ipoint_in=1
      if (mype/=ipe_now) xf(1,:)=zero
    else
      if (mype==0) then
        call mpistop('Field tracing error: given point is not in simulation box!')
      endif
    endif
 
 
    ! other points in field line
    do while(stopt .eqv. .false.)
 
      if (mype==ipe_now) then
        igrid=igrid_now
        ! looking for points in one pe
        call find_points_in_pe(igrid,ipoint_in,xf,wp,wl,dl,nump,nwp,nwl,forward,ftype,tcondi,statusf)
        status_bcast(1:4+ndim)=statusf(1:4+ndim)
        status_bcast(4+ndim+1:4+ndim+nwl)=wl(2:1+nwl)
      endif
      ! comunication
      call mpi_bcast(status_bcast,4+ndim+nwl,mpi_double_precision,ipe_now,icomm,ierrmpi)
      statusf(1:4+ndim)=status_bcast(1:4+ndim)
      wl(2:1+nwl)=status_bcast(4+ndim+1:4+ndim+nwl)
 
      ! prepare for next step
      ipoint_out=int(statusf(1))
      ipe_next=int(statusf(2))
      igrid_next=int(statusf(3))
      if (int(statusf(4))==1) then
        stopt=.true.
        wl(1)=ipoint_out-1
      endif
      if (mype==ipe_next) then
        do j=1,ndim
          xf(ipoint_out,j)=statusf(4+j)
        enddo
      else
        xf(ipoint_out,:)=zero
      endif
 
      ! pe and grid of next point
      ipe_now=ipe_next
      igrid_now=igrid_next
      ipoint_in=ipoint_out
    enddo
 
    if (tcondi/='TRAC') then
      numrt=int(wl(1))
      allocate(data_send(numrt,ndim+nwp),data_recv(numrt,ndim+nwp))
      data_send(:,:)=zero
      data_recv(:,:)=zero
      data_send(1:numrt,1:ndim)=xf(1:numrt,1:ndim)
      if (nwp>0) data_send(1:numrt,1+ndim:ndim+nwp)=wp(1:numrt,1:nwp)
      call mpi_allreduce(data_send,data_recv,numrt*(ndim+nwp),mpi_double_precision,&
                         mpi_sum,icomm,ierrmpi)
      xf(1:numrt,1:ndim)=data_recv(1:numrt,1:ndim)
      if (nwp>0) wp(1:numrt,1:nwp)=data_recv(1:numrt,1+ndim:ndim+nwp)
      deallocate(data_send,data_recv)
    endif
 
  end subroutine trace_field_single
 
  subroutine find_points_in_pe(igrid,ipoint_in,xf,wP,wL,dL,numP,nwP,nwL,forward,ftype,tcondi,statusF)
 
    integer, intent(inout) :: igrid
    integer, intent(in) :: ipoint_in,numP,nwP,nwL
    double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
    double precision, intent(in) :: dL
    logical, intent(in) :: forward
    character(len=std_len), intent(in) :: ftype,tcondi
    double precision, intent(inout) :: statusF(4+ndim)
 
    integer :: ipe_next,igrid_next,ip_in,ip_out,j,indomain
    logical :: newpe,stopT
    double precision :: xfout(ndim)
 
    ip_in=ipoint_in
    newpe=.false.
 
    do while(newpe .eqv. .false.)
      ! looking for points in given grid    
      call find_points_interp(igrid,ip_in,ip_out,xf,wp,wl,nump,nwp,nwl,dl,forward,ftype,tcondi)
      ip_in=ip_out
 
      ! when next point is out of given grid, find next grid  
      indomain=0
      {if (xf(ip_out,^db)>=xprobmin^db .and. xf(ip_out,^db)<=xprobmax^db) indomain=indomain+1\}      
      if (ip_out<nump .and. indomain==ndim) then
        if (tcondi/='TRAC') then
          stopt=.false.
          xfout=xf(ip_out,:)
          call find_next_grid(igrid,igrid_next,ipe_next,xfout,newpe,stopt)
        else
          if (xf(ip_out,ndim)>phys_trac_mask) then
            newpe=.true.
            stopt=.true.
          else
            stopt=.false.
            xfout=xf(ip_out,:)
            call find_next_grid(igrid,igrid_next,ipe_next,xfout,newpe,stopt)
          endif
        endif
      else
        newpe=.true.
        stopt=.true.
      endif
 
      if (newpe) then
        statusf(1)=ip_out
        statusf(2)=ipe_next
        statusf(3)=igrid_next
        statusf(4)=0
        if (stopt) statusf(4)=1
        do j=1,ndim
          statusf(4+j)=xf(ip_out,j)
        enddo
      endif
 
      if (newpe .eqv. .false.) igrid=igrid_next
    enddo
 
  end subroutine find_points_in_pe
 
  subroutine find_next_grid(igrid,igrid_next,ipe_next,xf1,newpe,stopT)
    ! check the grid and pe of next point
    use mod_usr_methods
    use mod_global_parameters
    use mod_forest
 
    integer, intent(inout) :: igrid,igrid_next,ipe_next
    double precision, intent(in) :: xf1(ndim)
    logical, intent(inout) :: newpe,stopT
 
    double precision :: dxb^D,xb^L,xbmid^D
    integer :: idn^D,my_neighbor_type,inblock
    integer :: ic^D,inc^D,ipe_neighbor,igrid_neighbor
    double precision :: xbn^L
 
    ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
    ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
    inblock=0
 
    ! direction of next grid
    idn^d=0\
    {if (xf1(^d)<=xbmin^d) idn^d=-1\}
    {if (xf1(^d)>=xbmax^d) idn^d=1\}
    my_neighbor_type=neighbor_type(idn^d,igrid)
    igrid_neighbor=neighbor(1,idn^d,igrid)
    ipe_neighbor=neighbor(2,idn^d,igrid)
 
    ! ipe and igrid of next grid
    select case(my_neighbor_type)
    case (neighbor_boundary)
      ! next point is not in simulation box
      newpe=.true.
      stopt=.true.
 
    case(neighbor_coarse)
      ! neighbor grid has lower refinement level      
      igrid_next=igrid_neighbor
      ipe_next=ipe_neighbor
      if (mype==ipe_neighbor) then
        newpe=.false.
      else
        newpe=.true.
      endif
 
    case(neighbor_sibling)
      ! neighbor grid has lower refinement level 
      igrid_next=igrid_neighbor
      ipe_next=ipe_neighbor
      if (mype==ipe_neighbor) then
        newpe=.false.
      else
        newpe=.true.
      endif
 
    case(neighbor_fine)
      ! neighbor grid has higher refinement level 
      {xbmid^d=(xbmin^d+xbmax^d)/2.d0\}
      ^d&inc^d=1\
      {if (xf1(^d)<=xbmin^d) inc^d=0\}
      {if (xf1(^d)>xbmin^d .and. xf1(^d)<=xbmid^d) inc^d=1\}
      {if (xf1(^d)>xbmid^d .and. xf1(^d)<xbmax^d) inc^d=2\}
      {if (xf1(^d)>=xbmax^d) inc^d=3\}
      ipe_next=neighbor_child(2,inc^d,igrid)
      igrid_next=neighbor_child(1,inc^d,igrid)
      if (mype==ipe_next) then
        newpe=.false.
      else
        newpe=.true.
      endif
    end select
 
  end subroutine find_next_grid
 
  subroutine find_points_interp(igrid,ip_in,ip_out,xf,wP,wL,numP,nwP,nwL,dL,forward,ftype,tcondi)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in) :: igrid,ip_in,numP,nwP,nwL
    integer, intent(inout) :: ip_out
    double precision, intent(inout) :: xf(numP,ndim),wP(numP,nwP),wL(1+nwL)
    double precision, intent(in) :: dL
    logical, intent(in) :: forward
    character(len=std_len), intent(in) :: ftype,tcondi
 
    double precision :: dxb^D,xb^L
    integer          :: ip,inblock,ixI^L,ixO^L,j
    double precision :: field(ixg^T,ndir)
    double precision :: xs1(ndim),xs2(ndim),K1(ndim),K2(ndim)
    double precision :: xfpre(ndim),xfnow(ndim),xfnext(ndim)
    double precision :: Tpre,Tnow,Tnext,dTds,Lt,Lr,ds,T_bott,trac_delta
 
    ixi^l=ixg^ll;
    ixo^l=ixm^ll;
    ^d&dxb^d=rnode(rpdx^d_,igrid);
    ^d&xbmin^d=rnode(rpxmin^d_,igrid);
    ^d&xbmax^d=rnode(rpxmax^d_,igrid);
 
    if (tcondi/='TRAC') then
      ds=dl
    else
      ds=dxb^nd
      t_bott=2.d4/unit_temperature
      trac_delta=0.25d0
    endif
 
    ! main loop
    mainloop: do ip=ip_in,nump-1
 
      ! integrate magnetic field with Runge-Kutta method
      xs1(:)=xf(ip,:)
      call get_k(xs1,igrid,k1,ixi^l,dxb^d,ftype)
      if (forward) then
        xs2(:)=xf(ip,:)+ds*k1(:)
      else
        xs2(:)=xf(ip,:)-ds*k1(:)
      endif
      call get_k(xs2,igrid,k2,ixi^l,dxb^d,ftype)
      if (forward) then
        xf(ip+1,:)=xf(ip,:)+ds*(0.5*k1(:)+0.5*k2(:))
      else
        xf(ip+1,:)=xf(ip,:)-ds*(0.5*k1(:)+0.5*k2(:))
      endif
      ip_out=ip+1
 
      ! get local values for variable via interpolation
      if (tcondi/='TRAC') then
        if (associated(usr_set_field_w)) then 
          call usr_set_field_w(igrid,ip,xf,wp,wl,nump,nwp,nwl,dl,forward,ftype,tcondi)
        endif
      else  ! get TRAC Tcoff
        wp(ip,1)=mype
        wp(ip,2)=igrid
        if (ip==ip_in) then
          if (forward) then
            xfpre(:)=xf(ip,:)-ds*k1(:)
          else
            xfpre(:)=xf(ip,:)+ds*k1(:)
          endif
          xfnow(:)=xf(ip,:)
          xfnext(:)=xf(ip+1,:)
          call get_t_loc_trac(igrid,xfpre,tpre,ixi^l,dxb^d)
          call get_t_loc_trac(igrid,xfnow,tnow,ixi^l,dxb^d)
          call get_t_loc_trac(igrid,xfnext,tnext,ixi^l,dxb^d)
        else
          xfpre=xf(ip-1,:)
          xfnow(:)=xf(ip,:)
          xfnext(:)=xf(ip+1,:)
          tpre=tnow
          tnow=tnext
          call get_t_loc_trac(igrid,xfnext,tnext,ixi^l,dxb^d)
        endif
        dtds=abs(tnext-tpre)/(2*ds)
        if (ip==1) then
          lt=0.d0
          wl(2)=t_bott   ! current Tcofl
          wl(3)=tnow     ! current Tlmax
        else
          lt=0.d0
          if (dtds>0.d0) then
            lt=tnow/dtds
            lr=ds
            ! renew cutoff temperature
            if(lr>trac_delta*lt) then
              if (tnow>wl(2)) wl(2)=tnow
            endif
          endif
          if (tnow>wl(3)) wl(3)=tnow
        endif
      endif
 
      ! exit loop if next point is not in this block
      inblock=0
      {if (xf(ip+1,^db)>=xbmin^db .and. xf(ip+1,^db)<xbmax^db) inblock=inblock+1\}
      if (tcondi=='TRAC' .and. xf(ip+1,ndim)>phys_trac_mask) inblock=0
      if (inblock/=ndim) exit mainloop
 
    enddo mainloop
 
  end subroutine find_points_interp
 
  subroutine get_k(xfn,igrid,K,ixI^L,dxb^D,ftype)
    use mod_usr_methods
 
    integer :: ixI^L,igrid
    double precision :: dxb^D
    double precision :: xfn(ndim),K(ndim)
    character(len=std_len) :: ftype
 
    integer          :: ixb^D,ix^D,ixbl^D,ixO^L,j
    double precision :: xd^D
    double precision :: field(0:1^D&,ndir),Fx(ndim),factor(0:1^D&)
    double precision :: vector(ixI^S,1:ndir)
    double precision :: Ftotal
 
    ^d&ixbl^d=floor((xfn(^d)-ps(igrid)%x(iximin^dd,^d))/dxb^d)+iximin^d;
    ^d&xd^d=(xfn(^d)-ps(igrid)%x(ixbl^dd,^d))/dxb^d;
    ^d&ixomin^d=ixbl^d;
    ^d&ixomax^d=ixbl^d+1;
    vector=zero
 
    field=zero
    if (ftype=='Bfield') then
      if (b0field) then
        vector(ixo^s,1:ndir)=ps(igrid)%w(ixo^s,iw_mag(1:ndir))+ps(igrid)%B0(ixo^s,1:ndir,0)
      else
        vector(ixo^s,1:ndir)=ps(igrid)%w(ixo^s,iw_mag(1:ndir))
      endif
    else if (ftype=='Vfield') then
      call phys_get_v(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,vector)
    endif
    {do ix^d=0,1\}
      factor(ix^d)={abs(1-ix^d-xd^d)*}
      field(ix^d,1:ndir)=vector(ixbl^d+ix^d,1:ndir)
    {enddo\}    
 
    fx=0.d0
    {do ix^db=0,1\}
      do j=1,ndim
        fx(j)=fx(j)+field(ix^d,j)*factor(ix^d)
      enddo
    {enddo\}
 
    if (ftype=='Bfield' .or. ftype=='Vfield') then
      fx=0.d0
      {do ix^db=0,1\}
        do j=1,ndim
          fx(j)=fx(j)+field(ix^d,j)*factor(ix^d)
        enddo
      {enddo\}
    else if (associated(usr_set_field)) then
      call usr_set_field(xfn,igrid,fx,ftype)
    else
      call mpistop('Field tracing error: wrong field type!')
    endif
 
    ftotal=zero
    do j=1,ndim
      ftotal=ftotal+(fx(j))**2
    enddo
    ftotal=dsqrt(ftotal)
 
    if (ftotal==0.d0) then
      k=0.d0
      k(1)=1.d0
    else
      k(1:ndim)=fx(1:ndim)/ftotal
    endif
 
  end subroutine get_k
 
  subroutine get_t_loc_trac(igrid,xloc,Tloc,ixI^L,dxb^D)
    ! grid T has been calculated and stored in wextra(ixI^S,Tcoff_)
    ! see mod_trac
    integer, intent(in) :: igrid,ixI^L
    double precision, intent(inout) :: xloc(ndim)
    double precision, intent(inout) :: Tloc
    double precision, intent(in) :: dxb^D
 
    integer          :: ixb^D,ix^D,ixbl^D,j,ixO^L
    double precision :: xd^D
    double precision :: factor(0:1^D&),Tnear(0:1^D&)
 
    ^d&ixbl^d=floor((xloc(^d)-ps(igrid)%x(iximin^dd,^d))/dxb^d)+iximin^d;
    ^d&xd^d=(xloc(^d)-ps(igrid)%x(ixbl^dd,^d))/dxb^d;
    ^d&ixomin^d=ixbl^d;
    ^d&ixomax^d=ixomin^d+1;
 
    {do ix^d=0,1\}
      factor(ix^d)={abs(1-ix^d-xd^d)*}
      tnear(ix^d)=ps(igrid)%wextra(ixbl^d+ix^d,iw_tcoff)
    {enddo\}
 
    tloc=0.d0
    ! interpolation
    {do ix^db=0,1\}
      tloc=tloc+tnear(ix^d)*factor(ix^d)
    {enddo\}
 
  end subroutine get_t_loc_trac
 
end module mod_trace_field