!---------------------------------------------------------------------------
SUBROUTINE ttheta_apply(nchir,ndim_tc,phi3,y3,ifirst)
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE APH_Module
USE Aph_Ham_Module
USE Cheby_Module
IMPLICIT NONE
!==============================================================================
! Apply T_theta:
!==============================================================================
! I N P U T / O U T P U T
INTEGER,INTENT(IN) :: nchir, ndim_tc,ifirst
REAL(dp),INTENT(IN)    :: phi3(nchir,nrho,ntheta)
REAL(dp),INTENT(INOUT) :: y3(nchir,nrho,ntheta)
!==============================================================================
! I N T E R N A L S
CHARACTER(len=100) :: filename
INTEGER :: irho, jrho, itheta, ichi, ii, iitheta, jjtheta
INTEGER :: prev_num, k
REAL(dp) :: term
!==============================================================================

 t_theta_clip=t_theta

 IF (clip_theta) THEN

   IF (first_theta) THEN

     prev_num=ntheta
     tempeig2=0.d0
     tempvcut2=0.d0
     num_eig2=0
     compr2=0.d0
     compr_mat2=0.d0

     DO ii=1,ntheta
       DO iitheta=1,ntheta
         DO jjtheta=1,ntheta
           compr2(jjtheta,iitheta)=tht_evec(jjtheta,ii)*tht_evec(iitheta,ii)
         ENDDO
       ENDDO
       compr_mat2(:,:,ii)=compr2
     ENDDO

! Determine how many eigenvales to replace at each rho and theta point

     filename='eignum_theta.txt'
     OPEN(110,file=TRIM(outdir)//TRIM(filename))
     DO irho=1,nrho
       k=0
       DO itheta=1,ntheta
         term=tht_eig(itheta)*rhom2(irho)
         IF (term.le.vcut) k=k+1
       ENDDO
       num_eig2(irho)=k
       write(110,*) irho, k
     ENDDO
     CLOSE(110)

     OPEN (unit=51,file=TRIM(BinOutdir)//'t_theta_clip.bin',form='unformatted')

    DO irho=1,nrho ! rho_loop

      jrho=nrho-irho+1

      IF (num_eig2(jrho).lt.prev_num) THEN
        DO ii=num_eig2(jrho+1),prev_num !ii loop

               tempeig2 =tempeig2 +(compr_mat2(:,:,ii)*tht_eig(ii))
               tempvcut2=tempvcut2+(compr_mat2(:,:,ii)*vcut)

        ENDDO ! end ii loop

        t_theta_clip=t_theta+(tempvcut2/rhom2(jrho))-tempeig2

      ELSEIF (num_eig2(jrho).eq.prev_num) THEN

        t_theta_clip=t_theta+(tempvcut2/rhom2(jrho))-tempeig2

      ENDIF

      IF (num_eig2(jrho).ne.ntheta) write(51) t_theta_clip

      y3(:,jrho,:)=Matmul(phi3(:,jrho,:),Transpose(t_theta_clip))*rhom2(jrho)

      !Call DGEMM('n','t',nchir,ntheta,ntheta,rhom2(irho),phi3(:,irho,:),
      !            nchir,t_theta,ntheta,0.d0,y3(:,irho,:), nchir)

      prev_num=num_eig2(jrho)

    ENDDO ! end rho loop

    first_theta=.false.

   ELSE

     OPEN (unit=51,file=TRIM(BinOutdir)//'t_theta_clip.bin',form='unformatted')

     DO irho=1,nrho ! rho_loop

!       t_theta_clip=t_theta

       jrho=nrho-irho+1

       IF (num_eig2(jrho).ne.ntheta) read(51) t_theta_clip

       y3(:,jrho,:)=Matmul(phi3(:,jrho,:),Transpose(t_theta_clip))*rhom2(jrho)

       !Call DGEMM('n','t',nchir,ntheta,ntheta,rhom2(irho),phi3(:,irho,:),
       !           nchir,t_theta,ntheta,0.d0,y3(:,irho,:),nchir)

     ENDDO ! end rho loop

     CLOSE(51)       

   ENDIF

ELSE

  ! Apply T_theta:
  DO irho=1,nrho !this is used elsewhere
      y3(:,irho,:)=Matmul(phi3(:,irho,:),Transpose(t_theta))*rhom2(irho)
  ENDDO

ENDIF

RETURN
END SUBROUTINE ttheta_apply