calc_grid Subroutine

  subroutine calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,&
     normconv,ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
     ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,first)
    ! this subroutine computes both corner as well as cell-centered values
    ! it handles how we do the center to corner averaging, as well as
    ! whether we switch to cartesian or want primitive or conservative output,
    ! handling the addition of B0 in B0+B1 cases, ...
    !
    ! the normconv is passed on to usr_aux_output for extending with
    ! possible normalization values for the nw+1:nw+nwauxio entries

    use mod_usr_methods, only: usr_aux_output
    use mod_global_parameters
    use mod_limiter
    use mod_physics, only: phys_energy, physics_type, phys_to_primitive
    use mod_comm_lib, only: mpistop

    integer, intent(in) :: qunit, igrid
    double precision, intent(in), dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
       ixMlo3-1:ixMhi3,ndim) :: xC
    double precision, intent(in), dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,&
       ixMlo3:ixMhi3,ndim)   :: xCC
    integer :: ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
       ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3
    logical, intent(in) :: first

    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
       ixMlo3-1:ixMhi3,ndim) :: xC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
       ndim)   :: xCC_TMP
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
       ixMlo3-1:ixMhi3,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
       nw+nwauxio)     :: wCC_TMP
    double precision,dimension(0:nw+nwauxio),intent(out)       :: normconv

    double precision :: ldw(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3),&
        dwC(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3)
    double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
       ixMlo3-1:ixMhi3,nw+nwauxio)   :: wC
    double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
       nw+nwauxio)     :: wCC
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw+nwauxio)   :: w
    integer :: nxCC1,nxCC2,nxCC3,idims,jxCmin1,jxCmin2,jxCmin3,jxCmax1,jxCmax2,&
       jxCmax3,iwe
    integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3, ix, iw, level, idir
    logical, save :: subfirst=.true.

    ! initialize w
    w=0.d0

    ixCmin1=ixMlo1-1;ixCmin2=ixMlo2-1;ixCmin3=ixMlo3-1; ixCmax1=ixMhi1
    ixCmax2=ixMhi2;ixCmax3=ixMhi3; !Corner indices
    ixCCmin1=ixMlo1;ixCCmin2=ixMlo2;ixCCmin3=ixMlo3; ixCCmax1=ixMhi1
    ixCCmax2=ixMhi2;ixCCmax3=ixMhi3; !Center indices

    nx1=ixMhi1-ixMlo1+1;nx2=ixMhi2-ixMlo2+1;nx3=ixMhi3-ixMlo3+1;
    level=node(plevel_,igrid)

    normconv(0) = length_convert_factor
    normconv(1:nw) = w_convert_factor
    block=>ps(igrid)
    w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw)=ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:nw)

    if (nwextra>0) then
     ! here we actually fill the ghost layers for the nwextra variables using
     ! continuous extrapolation (as these values do not exist normally in ghost cells)
     do idims=1,ndim
      select case(idims)
       case(1)
         jxCmin1=ixGhi1+1-nghostcells;jxCmin2=ixGlo2;jxCmin3=ixGlo3;
         jxCmax1=ixGhi1;jxCmax2=ixGhi2;jxCmax3=ixGhi3;
         do ix1=jxCmin1,jxCmax1
             w(ix1,jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw) = w(jxCmin1-1,jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw)
         end do
         jxCmin1=ixGlo1;jxCmin2=ixGlo2;jxCmin3=ixGlo3;
         jxCmax1=ixGlo1-1+nghostcells;jxCmax2=ixGhi2;jxCmax3=ixGhi3;
         do ix1=jxCmin1,jxCmax1
             w(ix1,jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw) = w(jxCmax1+1,jxCmin2:jxCmax2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw)
         end do 
       case(2)
         jxCmin1=ixGlo1;jxCmin2=ixGhi2+1-nghostcells;jxCmin3=ixGlo3;
         jxCmax1=ixGhi1;jxCmax2=ixGhi2;jxCmax3=ixGhi3;
         do ix2=jxCmin2,jxCmax2
             w(jxCmin1:jxCmax1,ix2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw) = w(jxCmin1:jxCmax1,jxCmin2-1,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw)
         end do
         jxCmin1=ixGlo1;jxCmin2=ixGlo2;jxCmin3=ixGlo3;
         jxCmax1=ixGhi1;jxCmax2=ixGlo2-1+nghostcells;jxCmax3=ixGhi3;
         do ix2=jxCmin2,jxCmax2
             w(jxCmin1:jxCmax1,ix2,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw) = w(jxCmin1:jxCmax1,jxCmax2+1,jxCmin3:jxCmax3,&
                nw-nwextra+1:nw)
         end do 
       case(3)
         jxCmin1=ixGlo1;jxCmin2=ixGlo2;jxCmin3=ixGhi3+1-nghostcells;
         jxCmax1=ixGhi1;jxCmax2=ixGhi2;jxCmax3=ixGhi3;
         do ix3=jxCmin3,jxCmax3
             w(jxCmin1:jxCmax1,jxCmin2:jxCmax2,ix3,&
                nw-nwextra+1:nw) = w(jxCmin1:jxCmax1,jxCmin2:jxCmax2,jxCmin3-1,&
                nw-nwextra+1:nw)
         end do
         jxCmin1=ixGlo1;jxCmin2=ixGlo2;jxCmin3=ixGlo3;
         jxCmax1=ixGhi1;jxCmax2=ixGhi2;jxCmax3=ixGlo3-1+nghostcells;
         do ix3=jxCmin3,jxCmax3
             w(jxCmin1:jxCmax1,jxCmin2:jxCmax2,ix3,&
                nw-nwextra+1:nw) = w(jxCmin1:jxCmax1,jxCmin2:jxCmax2,jxCmax3+1,&
                nw-nwextra+1:nw)
         end do 
      end select
     end do
    end if

    ! next lines needed when usr_aux_output uses gradients
    ! and later on when dwlimiter2 is used
    dxlevel(1)=rnode(rpdx1_,igrid);dxlevel(2)=rnode(rpdx2_,igrid)
    dxlevel(3)=rnode(rpdx3_,igrid);
    if(nwauxio>0)then
      ! auxiliary io variables can be computed and added by user
      ! next few lines ensure correct usage of routines like divvector etc
      ! default (no) normalization for auxiliary variables
      normconv(nw+1:nw+nwauxio)=one

      if (.not. associated(usr_aux_output)) then
         call mpistop("usr_aux_output not defined")
      else
         call usr_aux_output(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
            ixMlo1-1,ixMlo2-1,ixMlo3-1,ixMhi1+1,ixMhi2+1,ixMhi3+1,w,&
            ps(igrid)%x,normconv)
      end if
    endif

    ! In case primitives to be saved: use primitive subroutine
    !  extra layer around mesh only needed when storing corner values and averaging
    if(saveprim.and.first) call phys_to_primitive(ixGlo1,ixGlo2,ixGlo3,ixGhi1,&
       ixGhi2,ixGhi3,ixMlo1-1,ixMlo2-1,ixMlo3-1,ixMhi1+1,ixMhi2+1,ixMhi3+1,&
       w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:nw),ps(igrid)%x)

    if(B0field) then
    ! B0+B1 split handled here
      if(.not.saveprim.and.phys_energy) then
        w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,iw_e)=w(ixGlo1:ixGhi1,&
           ixGlo2:ixGhi2,ixGlo3:ixGhi3,iw_e)+&
           0.5d0*sum(ps(igrid)%B0(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,:,&
           0)**2,dim=ndim+1) + sum(w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
           iw_mag(:))*ps(igrid)%B0(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,:,&
           0),dim=ndim+1)
      end if
      w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,iw_mag(:))=w(ixGlo1:ixGhi1,&
         ixGlo2:ixGhi2,ixGlo3:ixGhi3,iw_mag(:))+ps(igrid)%B0(ixGlo1:ixGhi1,&
         ixGlo2:ixGhi2,ixGlo3:ixGhi3,:,0)
    end if
    ! compute the cell-center values for w first
    ! cell center values obtained from mere copy
    wCC(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,ixCCmin3:ixCCmax3,&
       :)=w(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,ixCCmin3:ixCCmax3,:)

    ! compute the corner values for w now by averaging

    if(slab_uniform) then
       ! for slab symmetry: no geometrical info required
       do iw=1,nw+nwauxio
         do ix3=ixCmin3,ixCmax3
         do ix2=ixCmin2,ixCmax2
         do ix1=ixCmin1,ixCmax1
            wC(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3:ix3+1,&
               iw))/dble(2**ndim)
         end do
         end do
         end do
       end do
    else
       do iw=1,nw+nwauxio
         do ix3=ixCmin3,ixCmax3
         do ix2=ixCmin2,ixCmax2
         do ix1=ixCmin1,ixCmax1
           wC(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3:ix3+1,&
              iw)*ps(igrid)%dvolume(ix1:ix1+1,ix2:ix2+1,&
              ix3:ix3+1)) /sum(ps(igrid)%dvolume(ix1:ix1+1,ix2:ix2+1,&
              ix3:ix3+1))
         end do
         end do
         end do
       end do
    endif

    if(nocartesian) then
      ! keep the coordinate and vector components
      xC_TMP(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
         1:ndim)          = xC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
         1:ndim)
      wC_TMP(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
         1:nw+nwauxio)    = wC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,ixCmin3:ixCmax3,&
         1:nw+nwauxio)
      xCC_TMP(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,ixCCmin3:ixCCmax3,&
         1:ndim)        = xCC(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,&
         ixCCmin3:ixCCmax3,1:ndim)
      wCC_TMP(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,ixCCmin3:ixCCmax3,&
         1:nw+nwauxio)  = wCC(ixCCmin1:ixCCmax1,ixCCmin2:ixCCmax2,&
         ixCCmin3:ixCCmax3,1:nw+nwauxio)
    else
      ! do all conversions to cartesian coordinates and vector components
      ! start for the corner values
      call to_cartesian(xC_TMP,wC_TMP,ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,&
         ixCmax3,xC,wC)
      ! then cell center values
      call to_cartesian(xCC_TMP,wCC_TMP,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
         ixCCmax2,ixCCmax3,xCC,wCC)
    endif

    ! Warning: differentiate between idl/idlCC/tecplot/tecplotCC/vtu(B)/vtu(B)CC
    if(nwaux>0 .and. mype==0 .and. first.and.subfirst) then
      ! when corner values are computed and auxiliaries present: warn!
      if(convert_type=='idl'.or.convert_type=='tecplot' &
         .or.convert_type=='vtu'.or.convert_type=='vtuB') write(*,&
         *) 'Warning: also averaged auxiliaries within calc_grid'
      subfirst=.false.
    endif

  end subroutine calc_grid