Curl_Vector3D_t Subroutine

  subroutine Curl_Vector3D_t(this,curlf)
    implicit none
    class(Vector3D_t),intent(in) :: this
    real(prec),intent(out) :: curlf(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:3)
    ! Local
    integer    :: i,j,k,ii,idir,iel,ivar
    real(prec) :: dfds1,dfds2,dfds3
    real(prec) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:3,1:3)

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar, &
                  idir=1:3)

      dfds1 = 0.0_prec
      dfds2 = 0.0_prec
      dfds3 = 0.0_prec
      do ii = 1,this%N+1
        dfds1 = dfds1+this%interp%dMatrix(ii,i)*this%interior(ii,j,k,iel,ivar,idir)
        dfds2 = dfds2+this%interp%dMatrix(ii,j)*this%interior(i,ii,k,iel,ivar,idir)
        dfds3 = dfds3+this%interp%dMatrix(ii,k)*this%interior(i,j,ii,iel,ivar,idir)
      enddo
      df(i,j,k,iel,ivar,idir,1) = dfds1
      df(i,j,k,iel,ivar,idir,2) = dfds2
      df(i,j,k,iel,ivar,idir,3) = dfds3

    enddo

    do concurrent(i=1:this%N+1,j=1:this%N+1, &
                  k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar)

      curlf(i,j,k,iel,ivar,1) = (df(i,j,k,iel,ivar,3,2)-df(i,j,k,iel,ivar,2,3))
      curlf(i,j,k,iel,ivar,2) = (df(i,j,k,iel,ivar,1,3)-df(i,j,k,iel,ivar,3,1))
      curlf(i,j,k,iel,ivar,3) = (df(i,j,k,iel,ivar,2,1)-df(i,j,k,iel,ivar,1,2))

    enddo

  endsubroutine Curl_Vector3D_t