prolong_2nd Subroutine

  subroutine prolong_2nd(sCo,ixComin1,ixComin2,ixComin3,ixComax1,ixComax2,&
     ixComax3,sFi,dxCo1,dxCo2,dxCo3,xComin1,xComin2,xComin3,dxFi1,dxFi2,dxFi3,&
     xFimin1,xFimin2,xFimin3,igridCo,igridFi)
    use mod_physics, only: phys_to_conserved
    use mod_global_parameters
    use mod_geometry, only: sync_positions_host
    use mod_amr_fct, only: already_fine, prolong_2nd_stg

    integer, intent(in) :: ixComin1,ixComin2,ixComin3,ixComax1,ixComax2,&
       ixComax3, igridFi, igridCo
    double precision, intent(in) :: dxCo1,dxCo2,dxCo3, xComin1,xComin2,xComin3,&
        dxFi1,dxFi2,dxFi3, xFimin1,xFimin2,xFimin3
    type(state), intent(in)      :: sCo
    type(state), intent(inout)   :: sFi

    integer :: ixCo1,ixCo2,ixCo3, jxCo1,jxCo2,jxCo3, hxCo1,hxCo2,hxCo3, ixFi1,&
       ixFi2,ixFi3, ix1,ix2,ix3, idim, iw, ixCgmin1,ixCgmin2,ixCgmin3,ixCgmax1,&
       ixCgmax2,ixCgmax3, el
    double precision :: slopeL, slopeR, slopeC, signC, signR
    double precision :: slope(nw,ndim)
    double precision :: eta1,eta2,eta3
    logical :: fine_min1,fine_min2,fine_min3,fine_max1,fine_max2,fine_max3

    ixCgmin1=ixComin1;ixCgmin2=ixComin2;ixCgmin3=ixComin3;ixCgmax1=ixComax1
    ixCgmax2=ixComax2;ixCgmax3=ixComax3;

    
    do ixCo3 = ixCgmin3,ixCgmax3
       do ixCo2 = ixCgmin2,ixCgmax2
          do ixCo1 = ixCgmin1,ixCgmax1
             ixFi3=2*(ixCo3-ixComin3)+ixMlo3
             ixFi2=2*(ixCo2-ixComin2)+ixMlo2
             ixFi1=2*(ixCo1-ixComin1)+ixMlo1

             do idim=1,ndim
                hxCo1=ixCo1-kr(1,idim)
                hxCo2=ixCo2-kr(2,idim)
                hxCo3=ixCo3-kr(3,idim)
                jxCo1=ixCo1+kr(1,idim)
                jxCo2=ixCo2+kr(2,idim)
                jxCo3=ixCo3+kr(3,idim)

                do iw=1,nwgc  ! analytic extras past nwgc are set in alloc_node
                   slopeL=bg(1)%w(ixCo1,ixCo2,ixCo3,iw, igridCo) - bg(1)%w(hxCo1,hxCo2,hxCo3,iw, igridCo)
                   slopeR=bg(1)%w(jxCo1,jxCo2,jxCo3,iw, igridCo) - bg(1)%w(ixCo1,ixCo2,ixCo3,iw, igridCo)
                   slopeC=half*(slopeR+slopeL)

                   ! get limited slope
                   signR=sign(one,slopeR)
                   signC=sign(one,slopeC)
                   !select case(prolong_limiter)
                   !case(1)
                   !  ! unlimited
                   !  slope(iw,idim)=slopeC
                   !case(2)
                   !  ! minmod
                   !  slope(iw,idim)=signR*max(zero,min(dabs(slopeR), &
                   !                                    signR*slopeL))
                   !case(3)
                   !  ! woodward
                   !  slope(iw,idim)=two*signR*max(zero,min(dabs(slopeR), &
                   !                     signR*slopeL,signR*half*slopeC))
                   !case(4)
                   !  ! koren
                   !  slope(iw,idim)=signR*max(zero,min(two*signR*slopeL, &
                   !   (dabs(slopeR)+two*slopeL*signR)*third,two*dabs(slopeR)))
                   !case default
                   slope(iw,idim)=signC*max(zero,min(dabs(slopeC), signC*slopeL,&
                        signC*slopeR))
                   !end select
                end do
             end do
             ! cell-centered coordinates of coarse grid point
             !^D&xCo^D=xCo({ixCo^DD},^D)
             do ix3=ixFi3,ixFi3+1
                do ix2=ixFi2,ixFi2+1
                   do ix1=ixFi1,ixFi1+1
                      ! cell-centered coordinates of fine grid point
                      !^D&xFi^D=xFi({ix^DD},^D)
                      if(slab_uniform) then
                         ! normalized distance between fine/coarse cell center
                         ! in coarse cell: ranges from -0.5 to 0.5 in each direction
                         ! (origin is coarse cell center)
                         ! hence this is +1/4 or -1/4 on cartesian mesh
                         !eta^D=(xFi^D-xCo^D)*invdxCo^D;
                         eta1=0.5d0*(dble(ix1-ixFi1)-0.5d0)
                         eta2=0.5d0*(dble(ix2-ixFi2)-0.5d0)
                         eta3=0.5d0*(dble(ix3-ixFi3)-0.5d0);
                      else
                         ! forefactor is -0.5d0 when ix=ixFi and +0.5d0 for ixFi+1
                         eta1=(dble(ix1-ixFi1)-0.5d0)*(one-bgeo%dvolume(ix1,ix2,&
                              ix3,igridFi) /sum(bgeo%dvolume(ixFi1:ixFi1+1,ix2,ix3,igridFi)))
                         ! forefactor is -0.5d0 when ix=ixFi and +0.5d0 for ixFi+1
                         eta2=(dble(ix2-ixFi2)-0.5d0)*(one-bgeo%dvolume(ix1,ix2,&
                              ix3,igridFi) /sum(bgeo%dvolume(ix1,ixFi2:ixFi2+1,ix3,igridFi)))
                         ! forefactor is -0.5d0 when ix=ixFi and +0.5d0 for ixFi+1
                         eta3=(dble(ix3-ixFi3)-0.5d0)*(one-bgeo%dvolume(ix1,ix2,&
                              ix3,igridFi) /sum(bgeo%dvolume(ix1,ix2,ixFi3:ixFi3+1,igridFi)))
                      end if
                      bg(1)%w(ix1,ix2,ix3,1:nwgc, igridFi) = bg(1)%w(ixCo1,ixCo2,ixCo3,1:nwgc, igridCo) + (slope(1:nwgc,&
                           1)*eta1)+(slope(1:nwgc,2)*eta2)+(slope(1:nwgc,3)*eta3)
                   end do
                end do
             end do
          end do
       end do
    end do
    if(stagger_grid) then
       call already_fine(sFi,igridFi,fine_min1,fine_min2,fine_min3,fine_max1,&
            fine_max2,fine_max3)
       call prolong_2nd_stg(sCo,sFi,ixComin1,ixComin2,ixComin3,ixComax1,&
            ixComax2,ixComax3,ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,dxCo1,&
            dxCo2,dxCo3,xComin1,xComin2,xComin3,dxFi1,dxFi2,dxFi3,xFimin1,xFimin2,&
            xFimin3,.false.,fine_min1,fine_min2,fine_min3,fine_max1,fine_max2,&
            fine_max3)
    end if

    ! this reads sFi%x on the host
    if(prolongprimitive) call sync_positions_host(igridFi)
    if(prolongprimitive) call phys_to_conserved(ixGlo1,ixGlo2,ixGlo3,ixGhi1,&
         ixGhi2,ixGhi3,ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,sFi%w,sFi%x)

  end subroutine prolong_2nd