calculate_ittp Subroutine

  subroutine calculate_ittp(can_be_ttp_flag, is_ttp_value)
    use theta_grid, only: ntgrid
    implicit none
    logical, dimension(nlambda), intent(out) :: can_be_ttp_flag
    logical, dimension(-ntgrid:ntgrid, nlambda), intent(out) :: is_ttp_value
    integer, dimension(-ntgrid:ntgrid) :: ittp_indices
    integer :: ig, il
    ittp_indices = nlambda+1
    can_be_ttp_flag = .false.
    is_ttp_value = .false.
    if (.not. grid_has_trapped_particles()) return

    ! Note we exclude the possibility of totally trapped particles
    ! existing at either end of the theta grid.
    do ig = -ntgrid+1, ntgrid-1
       ! all pitch angles greater than or equal to ittp are totally trapped or forbidden
       do il = ng2+1, nlambda
          if (forbid(ig-1,il) .and. forbid(ig+1, il) &
               .and. .not. forbid(ig, il)) then
             ittp_indices(ig) = il

             ! Record the lowest il which satisfies the condition
             exit
          end if
       end do
    end do

    do il = ng2+1, nlambda
       can_be_ttp_flag(il) = any(il >= ittp_indices)
    end do

    ! Calculate is_ttp_value, indicating which pitch angles can reach
    ! and are totally trapped at this theta location.
    do il = ng2 + 1, nlambda
       do ig = -ntgrid, ntgrid
          if (il >= ittp_indices(ig) .and. .not. forbid(ig, il)) is_ttp_value(ig, il) = .true.
       end do
    end do

  end subroutine calculate_ittp