imk_initialize Subroutine

  subroutine imk_initialize(this, lev)
    class(pf_imk_t),   intent(inout) :: this
    class(pf_level_t), intent(inout) :: lev

    integer :: m, nnodes

    this%npieces = 1
    nnodes = lev%nnodes

    allocate(this%QdiffE(nnodes-1,nnodes), this%QtilE(nnodes-1,nnodes))
    allocate(this%dtsdc(nnodes-1))
    allocate(this%tsdc(nnodes))


    this%dtsdc = lev%nodes(2:nnodes) - lev%nodes(1:nnodes-1)
    this%bernoullis = 0.0_pfdp
    this%bernoullis(1 ) =       -1.0_pfdp / 2.0_pfdp
    this%bernoullis(2 ) =        1.0_pfdp / 6.0_pfdp
    this%bernoullis(4 ) =       -1.0_pfdp / 3.0_pfdp
    this%bernoullis(6 ) =        1.0_pfdp / 4.2e1_pfdp
    this%bernoullis(8 ) =       -1.0_pfdp / 3.0e1_pfdp
    this%bernoullis(10) =        5.0_pfdp / 6.6e1_pfdp
    this%bernoullis(12) =     -691.0_pfdp / 2.73e3_pfdp
    this%bernoullis(14) =        7.0_pfdp / 6.0_pfdp
    this%bernoullis(16) =    -3617.0_pfdp / 5.10e2_pfdp
    this%bernoullis(18) =    43867.0_pfdp / 7.98e2_pfdp
    this%bernoullis(20) =    -174611.0_pfdp/330.0_pfdp
    !>  Assign explicit approximate quadrature rule
    this%QtilE =  lev%sdcmats%qmatFE
    this%QdiffE = lev%sdcmats%qmat-this%QtilE

    !>  Make space for temporary variables
    call lev%ulevel%factory%create_array(this%A, nnodes, &
         lev%index,   lev%shape)

    do m = 1, nnodes
       call this%A(m)%setval(0.0_pfdp)
    end do

  end subroutine imk_initialize