text:
module ad_intrinsic

    use ad_kinds, only: rk
    use ad_types, only: tree_t, assignment(=)
    implicit none
    private

    public :: abs, max, min
    public :: sin, cos, tan
    public :: exp, sqrt, log, log10

    interface abs
        module procedure :: abs_t
    end interface abs

    interface max
        module procedure :: max_tt
        module procedure :: max_tr
        ! Ambiguous interfaces in generic interface 'max' for 'max_tr' at (1) and 'max_rt' here
        ! module procedure :: max_rt
    end interface max

    interface min
        module procedure :: min_tt
        module procedure :: min_tr
        ! module procedure :: min_rt
    end interface min

    interface sin
        module procedure :: sin_t
    end interface sin

    interface cos
        module procedure :: cos_t
    end interface cos

    interface tan
        module procedure :: tan_t
    end interface tan

    interface exp
        module procedure :: exp_t
    end interface exp

    interface sqrt
        module procedure :: sqrt_t
    end interface sqrt

    interface log
        module procedure :: log_t
    end interface log

    interface log10
        module procedure :: log10_t
    end interface log10

contains

    impure elemental function abs_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = abs(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = merge(-1.0_rk, 1.0_rk, t1%node%value < 0.0_rk)

    end function abs_t

    impure elemental function max_tt(t1, t2) result(t)
        type(tree_t), intent(in) :: t1, t2
        type(tree_t) :: t

        t = max(t1%node%value, t2%node%value)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk

    end function max_tt

    ! impure elemental function max_rt(r, t1) result(t)
    !     real(rk), intent(in) :: r
    !     type(tree_t), intent(in) :: t1
    !     type(tree_t) :: t

    !     allocate (t%node)
    !     t%node%value = max(r, t1%node%value)

    !     t%node%left => t1%node
    !     t%node%left_grad = 1.0_rk

    ! end function max_rt

    impure elemental function max_tr(t1, r) result(t)
        type(tree_t), intent(in) :: t1
        real(rk), intent(in) :: r
        type(tree_t) :: t

        t = max(t1%node%value, r)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk

    end function max_tr

    impure elemental function min_tt(t1, t2) result(t)
        type(tree_t), intent(in) :: t1, t2
        type(tree_t) :: t

        t = min(t1%node%value, t2%node%value)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk

    end function min_tt

    ! impure elemental function min_rt(r, t1) result(t)
    !     real(rk), intent(in) :: r
    !     type(tree_t), intent(in) :: t1
    !     type(tree_t) :: t

    !     allocate (t%node)
    !     t%node%value = min(r, t1%node%value)

    !     t%node%left => t1%node
    !     t%node%left_grad = 1.0_rk

    ! end function min_rt

    impure elemental function min_tr(t1, r) result(t)
        type(tree_t), intent(in) :: t1
        real(rk), intent(in) :: r
        type(tree_t) :: t

        t = min(t1%node%value, r)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk

    end function min_tr

    impure elemental function sin_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = sin(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = cos(t1%node%value)

    end function sin_t

    impure elemental function cos_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = cos(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = -sin(t1%node%value)

    end function cos_t

    impure elemental function tan_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = tan(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk/(cos(t1%node%value)*cos(t1%node%value))

    end function tan_t

    impure elemental function exp_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = exp(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = t%node%value

    end function exp_t

    impure elemental function sqrt_t(t1) result(t)
        use, intrinsic :: ieee_arithmetic, only: ieee_value, NAN => ieee_quiet_nan
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = sqrt(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = merge(0.5_rk/t%node%value, ieee_value(1.0_rk, NAN), t1%node%value >= 0.0_rk) !TODO: NAN

    end function sqrt_t

    impure elemental function log_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = log(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk/t1%node%value

    end function log_t

    impure elemental function log10_t(t1) result(t)
        type(tree_t), intent(in) :: t1
        type(tree_t) :: t

        t = log10(t1%node%value)

        t%node%left => t1%node
        t%node%left_grad = 1.0_rk/t1%node%value/log(10.0_rk)

    end function log10_t

end module ad_intrinsic

max_stars_repo_name: zoziha/Auto-Diff

max_issues_repo_name: zoziha/Auto-Diff

max_forks_repo_name: zoziha/Auto-Diff


text:
      subroutine qqb_w2jet_gvecx(p,n,in,msq,msqv_cs,msqvx)
      implicit none
c----Matrix element for W+2jet production
C----averaged over initial colours and spins
c    line 6 contracted with the vector n(mu)
C For nwz=+1
c     u(-p1)+dbar(-p2)--> g(p5)+ g(p6)+W^+(n(p3)+e^+(p4))
C For nwz=-1
c     d(-p1)+ubar(-p2)--> g(p5)+ g(p6)+W^-(e^-(p3)+nbar(p4))
c---It has been checked that this gives the right matrix element
c---squared when n is replaced by two physical polarizations,
c---but this routine contains one factor of 1/2 for identical gluons
c   in the final state.
      include 'constants.f'
      include 'masses.f'
      include 'qcdcouple.f'
      include 'ewcouple.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'ckm.f'
C ip is the label of the emitter
      integer j,k,in,i,n1,n2
      double precision msq(-nf:nf,-nf:nf),p(mxpart,4)
      double precision fac,prop,n(4),Vfac
      double precision p1p2(0:2,-1:1,-1:1),q1q2(0:2,-1:1,-1:1)
      double complex zab(mxpart,mxpart),zba(mxpart,mxpart)
      double precision msqv_cs(0:2,-nf:nf,-nf:nf)
      double precision msqvx(-nf:nf,-nf:nf,-nf:nf,-nf:nf)
      common/p1p2/p1p2
      common/q1q2/q1q2
!$omp threadprivate(/p1p2/,/q1q2/)
c--- note that we will use the first index of p1p2 to label
c--- the colour structure of the squared matrix element:
c---     0 --> -ninth*(qed piece)
c---     1 --> +(qcd ordering 1)
c---     2 --> +(qcd ordering 1)
c---     3 --> Total (0+1+2)
c--- In this (x) version, the array q1q2 is supposed to provide
c--- the same information as p1p2, but with the identities of the
c--- two outgoing particles switched (eg. gg->qqb rather than gg->qbq)

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      do i=0,2
        msqv_cs(i,j,k)=0d0
      enddo
      enddo
      enddo

      do i=0,2
      do j=-1,1
      do k=-1,1
      q1q2(i,j,k)=0d0
      p1p2(i,j,k)=0d0
      enddo
      enddo
      enddo

      prop=s(3,4)**2/((s(3,4)-wmass**2)**2+wmass**2*wwidth**2)
      fac=V*xn/eight*(gsq*gwsq)**2*prop
      call spinoru(6,p,za,zb)
C   zab=<i-|k|j-> zba=<i+|k|j+> where k is an arbitrary 4-vector
c---Arbitrary conventions of Bern, Dixon, Kosower, Weinzierl,
c---za(i,j)*zb(j,i)=s(i,j)
      call spinork(6,p,zab,zba,n)

      if (in .eq. 1) then
C--initial-initial
        call w2jetnx(2,5,3,4,6,1,p,n,za,zb,zab,zba)
        call storecsv_px(0,+1)
        call w2jetnx(5,2,3,4,6,1,p,n,za,zb,zab,zba)
        call storecsv_px(0,-1)
c--- g g -> qb q
        call w2jetnx(5,6,3,4,2,1,p,n,za,zb,zab,zba)
        call storecsv_px(0,0)
c--- g g -> q qb
        call w2jetnx(6,5,3,4,2,1,p,n,za,zb,zab,zba)
        call storecsv_qx(0,0)
        do i=0,2
          p1p2(i,0,+1)=aveqg*fac*p1p2(i,0,+1)
          p1p2(i,0,-1)=aveqg*fac*p1p2(i,0,-1)
          p1p2(i,0, 0)=avegg*fac*p1p2(i,0,0)
          q1q2(i,0, 0)=avegg*fac*q1q2(i,0,0)
        enddo
c      p1p2(0,+1)=aveqg*fac*w2jetn(2,6,3,4,5,1,p,n,za,zb,zab,zba)
c      p1p2(0,-1)=aveqg*fac*w2jetn(6,2,3,4,5,1,p,n,za,zb,zab,zba)
      elseif (in .eq. 2) then
        call w2jetnx(1,5,3,4,6,2,p,n,za,zb,zab,zba)
        call storecsv_px(+1,0)
        call w2jetnx(5,1,3,4,6,2,p,n,za,zb,zab,zba)
        call storecsv_px(-1,0)
c--- g g -> qb q
        call w2jetnx(5,6,3,4,1,2,p,n,za,zb,zab,zba)
        call storecsv_px(0,0)
c--- g g -> q qb
        call w2jetnx(6,5,3,4,1,2,p,n,za,zb,zab,zba)
        call storecsv_qx(0,0)
        do i=0,2
          p1p2(i,+1,0)=aveqg*fac*p1p2(i,+1,0)
          p1p2(i,-1,0)=aveqg*fac*p1p2(i,-1,0)
          p1p2(i,0, 0)=avegg*fac*p1p2(i,0,0)
          q1q2(i,0, 0)=avegg*fac*q1q2(i,0,0)
        enddo
c      p1p2(+1,0)=aveqg*fac*w2jetn(1,6,3,4,5,2,p,n,za,zb,zab,zba)
c      p1p2(-1,0)=aveqg*fac*w2jetn(6,1,3,4,5,2,p,n,za,zb,zab,zba)
      elseif (in .eq. 5) then
        call w2jetnx(1,2,3,4,6,5,p,n,za,zb,zab,zba)
        call storecsv_px(1,-1)
        call w2jetnx(2,1,3,4,6,5,p,n,za,zb,zab,zba)
        call storecsv_px(-1,1)
        call w2jetnx(1,6,3,4,2,5,p,n,za,zb,zab,zba)
        call storecsv_px(+1,0)
        call w2jetnx(6,1,3,4,2,5,p,n,za,zb,zab,zba)
        call storecsv_px(-1,0)
        call w2jetnx(2,6,3,4,1,5,p,n,za,zb,zab,zba)
        call storecsv_px(0,+1)
        call w2jetnx(6,2,3,4,1,5,p,n,za,zb,zab,zba)
        call storecsv_px(0,-1)
        do i=0,2
          p1p2(i,1,-1)=half*aveqq*fac*p1p2(i,1,-1)
          p1p2(i,-1,1)=half*aveqq*fac*p1p2(i,-1,1)
          p1p2(i,+1,0)=aveqg*fac*p1p2(i,+1,0)
          p1p2(i,-1,0)=aveqg*fac*p1p2(i,-1,0)
          p1p2(i,0,+1)=aveqg*fac*p1p2(i,0,+1)
          p1p2(i,0,-1)=aveqg*fac*p1p2(i,0,-1)
        enddo
c      p1p2(1,-1)=aveqq*fac*w2jetn(1,2,3,4,6,5,p,n,za,zb,zab,zba)
c      p1p2(-1,1)=aveqq*fac*w2jetn(2,1,3,4,6,5,p,n,za,zb,zab,zba)
c      p1p2(+1,0)=aveqg*fac*w2jetn(1,6,3,4,2,5,p,n,za,zb,zab,zba)
c      p1p2(-1,0)=aveqg*fac*w2jetn(6,1,3,4,2,5,p,n,za,zb,zab,zba)
c      p1p2(0,+1)=aveqg*fac*w2jetn(2,6,3,4,1,5,p,n,za,zb,zab,zba)
c      p1p2(0,-1)=aveqg*fac*w2jetn(6,2,3,4,1,5,p,n,za,zb,zab,zba)
      elseif (in .eq. 6) then
        call w2jetnx(1,2,3,4,5,6,p,n,za,zb,zab,zba)
        call storecsv_px(1,-1)
        call w2jetnx(2,1,3,4,5,6,p,n,za,zb,zab,zba)
        call storecsv_px(-1,1)
        call w2jetnx(1,5,3,4,2,6,p,n,za,zb,zab,zba)
        call storecsv_px(+1,0)
        call w2jetnx(5,1,3,4,2,6,p,n,za,zb,zab,zba)
        call storecsv_px(-1,0)
        call w2jetnx(2,5,3,4,1,6,p,n,za,zb,zab,zba)
        call storecsv_px(0,+1)
        call w2jetnx(5,2,3,4,1,6,p,n,za,zb,zab,zba)
        call storecsv_px(0,-1)
        do i=0,2
          p1p2(i,1,-1)=half*aveqq*fac*p1p2(i,1,-1)
          p1p2(i,-1,1)=half*aveqq*fac*p1p2(i,-1,1)
          p1p2(i,+1,0)=aveqg*fac*p1p2(i,+1,0)
          p1p2(i,-1,0)=aveqg*fac*p1p2(i,-1,0)
          p1p2(i,0,+1)=aveqg*fac*p1p2(i,0,+1)
          p1p2(i,0,-1)=aveqg*fac*p1p2(i,0,-1)
        enddo
c      p1p2(1,-1)=aveqq*fac*w2jetn(1,2,3,4,5,6,p,n,za,zb,zab,zba)
c      p1p2(-1,1)=aveqq*fac*w2jetn(2,1,3,4,5,6,p,n,za,zb,zab,zba)
c      p1p2(+1,0)=aveqg*fac*w2jetn(1,5,3,4,2,6,p,n,za,zb,zab,zba)
c      p1p2(-1,0)=aveqg*fac*w2jetn(5,1,3,4,2,6,p,n,za,zb,zab,zba)
c      p1p2(0,+1)=aveqg*fac*w2jetn(2,5,3,4,1,6,p,n,za,zb,zab,zba)
c      p1p2(0,-1)=aveqg*fac*w2jetn(5,2,3,4,1,6,p,n,za,zb,zab,zba)
      endif

      do j=-nf,nf
      do k=-nf,nf
      if     ((j .gt. 0) .and. (k .lt. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=Vsq(j,k)*p1p2(i,1,-1)
          enddo
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=Vsq(j,k)*p1p2(i,-1,1)
          enddo
      elseif ((j .gt. 0) .and. (k .eq. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=
     &  (Vsq(j,-1)+Vsq(j,-2)+Vsq(j,-3)+Vsq(j,-4)+Vsq(j,-5))*p1p2(i,+1,0)
          enddo
          do n1=1,nf
          msqvx(j,k,n1,0)=Vsq(j,-n1)*(
     .                         p1p2(0,+1,0)+p1p2(1,+1,0)+p1p2(2,+1,0))
          msqvx(j,k,0,n1)=msqvx(j,k,n1,0)
          enddo
      elseif ((j .lt. 0) .and. (k .eq. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=
     &  (Vsq(j,+1)+Vsq(j,+2)+Vsq(j,+3)+Vsq(j,+4)+Vsq(j,+5))*p1p2(i,-1,0)
          do n1=1,nf
          msqvx(j,k,-n1,0)=Vsq(j,n1)*(
     .                          p1p2(0,-1,0)+p1p2(1,-1,0)+p1p2(2,-1,0))
          msqvx(j,k,0,-n1)=msqvx(j,k,-n1,0)
          enddo
          enddo
      elseif ((j .eq. 0) .and. (k .gt. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=
     &  (Vsq(-1,k)+Vsq(-2,k)+Vsq(-3,k)+Vsq(-4,k)+Vsq(-5,k))*p1p2(i,0,+1)
          enddo
          do n1=1,nf
          msqvx(j,k,n1,0)=Vsq(-n1,k)*(
     .                          p1p2(0,0,+1)+p1p2(1,0,+1)+p1p2(2,0,+1))
          msqvx(j,k,0,n1)=msqvx(j,k,n1,0)
          enddo
      elseif ((j .eq. 0) .and. (k .lt. 0)) then
          do i=0,2
            msqv_cs(i,j,k)=
     &  (Vsq(+1,k)+Vsq(+2,k)+Vsq(+3,k)+Vsq(+4,k)+Vsq(+5,k))*p1p2(i,0,-1)
          enddo
          do n1=1,nf
          msqvx(j,k,-n1,0)=Vsq(n1,k)*(
     .                          p1p2(0,0,-1)+p1p2(1,0,-1)+p1p2(2,0,-1))
          msqvx(j,k,0,-n1)=msqvx(j,k,-n1,0)
          enddo
      elseif ((j .eq. 0) .and. (k .eq. 0)) then
          Vfac=0d0
          do n1=1,nf
            do n2=-nf,-1
              Vfac=Vfac+Vsq(n1,n2)
              msqvx(j,k,-n1,-n2)=Vsq(n1,n2)*(
     .                            p1p2(0,0,0)+p1p2(1,0,0)+p1p2(2,0,0))
              msqvx(j,k,-n2,-n1)=Vsq(n1,n2)*(
     .                            q1q2(0,0,0)+q1q2(1,0,0)+q1q2(2,0,0))
            enddo
          enddo
          do i=0,2
            msqv_cs(i,j,k)=Vfac*p1p2(i,0,0)
          enddo
      endif
      msq(j,k)=msqv_cs(0,j,k)+msqv_cs(1,j,k)+msqv_cs(2,j,k)
      enddo
      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
c***********************************************************************
c***********************************************************************
c dmm 20000330
c
c     genRhos subroutine
c
c     This subroutine calculates and loads values into the xrhoofn
c     tables, replacing the recursion and gamma-table based rhoOfN
c     calculations used previously.
c
c     The sums-of-states algorithm used, SBS in 'straightBS.f',
c     automatically convolves with any number of 1d-rotors, assuming
c     them classical after method of Astholz (see Gilbert and Smith).
c     IN this particular implementation, we can set all rotational
c     constants and symmetries to 1.0, since they will cancel (all
c     implementations are ratios).  Unit issues for the energies should
c     also cancel.
c
c     There is a lot of complicated machinery in here to deal with
c     noninteger degeneracies from the frequency fitting.  If we could
c     force integer frequencies, lots of this wouldn't be needed.
c
c     eMax and frequencies are expected in cm-1.  The grainsize should
c     be in cdnos.fh, its default set (10cm-1) or read in by cdgets.f
c
c     MODIFICATIONS
c
c     dmm 20000627
c     Adding 'wellscalcd' vector, which is initialized to
c     0.0 and stores the maximum temperature for which the xrhooofn
c     vector for a particular array was last calculated.  If new T is
c     greater than this T, recalc the well, otherwise, skip (should be
c     OK).  This should not affect the operation of the maxBSindex
c     vector, which lets us generate a warning when a call to xNOS
c     exceeds the energy levels to which we have calculated rho.
c
c***********************************************************************
c***********************************************************************

      subroutine genRhos(eMax, dE, temp)

      implicit real*8 (a-h, o-z)
      include 'cdparams.fh'
      include 'straightBS.fh'
      include 'cdnos.fh'
      include 'cdisprop.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdcontrl.fh'

      integer numfreq(mxWells, mxFreqs)
      integer iWell, Isize, imax1, ip1, iw1, nfrot, ndos
      real*8 VIB(maxmodes), BR(maxmodes), BDG(maxmodes)
      real*8 remaind(mxFreqs)
      real*8 Rcm
      real*8 Egrain1, Emax2, Viblo, Ebarmax
      integer maxnmax = 0

c     Set partition function flag.  See cdnos.fh.  Will be problematic
c     here if NOS's are pre-calc'd somewhere else

      pfflag = .false.
c

c     Boltz. constant cm-1/mol/K
      Rcm = 0.694995954

c     Set all BR, BDG to 1.0 as per discussion above
      do 40 i=1, maxmodes
         BR(i) = 20.0D0
c         BR(i) = 1.0D0
         BDG(i) = 1.0D0
 40   continue

c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** START                                     c
c----------------------------------------------------------------------c
c     Loop over wells; go one more if RRKM option is turned on to get
c     the RRKM TS sum of states.

      if(RRKM) then
         ndos = nWells + 1
      else
         ndos = nWells
      endif
c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** END                                       c
c----------------------------------------------------------------------c

      do 10 iWell = 1, ndos

c***********************************************************************
c     Begin If-statement block (goes all the way to '10 continue') which
c     is entered if the density of states for this well needs to be
c     calculated.
c***********************************************************************

c     dmm 20000908 this logic is screwed up and inpwonly actually isn't
c     available here anyway.

c     dmm 20010612 commented out below because I don't understand why '
c     inpwonly' logic is there ...
c$$$         if(temp.gt.wellscalcd(iWell) .or. (.not. inpwonly)) then
c$$$            wellscalcd(iWell) = temp
         if(temp.gt.wellscalcd(iWell)) then

c     set nfrot, number of free internal rotors, to 0
            nfrot = 0

c     2.  Find sum of degeneracies degensum
c     3.  Find the integer degeneracies for the BS
c     4.  Find the sum of the integer degeneracies ntotal
c     5.  Load the VIB vector for the BS algorithm

            sumdegen = 0.0
            ntotal = 0
            jcount = 1

            do 20 ifreq = 1, nFreqs(iWell)

               sumdegen = sumdegen + degen(iWell, ifreq)
               numfreq(iWell, ifreq) = nint(degen(iWell, ifreq))
               remaind(ifreq) = degen(iWell, ifreq) -
     $              dfloat(numfreq(iWell,ifreq))

               ntotal = ntotal + numfreq(iWell, ifreq)

c     Avoid loading the VIB vector if the frequencies have been read
c     from a file
               if(.not.FrRtRd(iWell)) then
                  do 15 j = 1, numfreq(iWell, ifreq)
                     VIB(jcount) = freq(iWell, ifreq)
                     jcount = jcount + 1
 15               continue
               endif
 20         continue

c***********************************************************************
c     Cleanup:  Check whether nint(sumdegen) is the same as summation of
c     frequences ntotal.  If not (if ntotal is less) then add that
c     frequency whose degeneracy is closest to the next integer number
c     .This complicated machinery is only necessary if we plan on using
c     noninteger degeneracies ...
c***********************************************************************

c     It shouldn't matter, but avoid this block if we have read the
c     frequencies from a file ...

            if(.not.FrRtRd(iWell)) then
               if( nint(sumdegen - 0.01).gt.ntotal) then
                  rmax = 0.0D0
c     Look for largest remainder
                  do 30 ifreq = 1, nFreqs(iWell)
                     if (remaind(ifreq).gt.rmax) then
                        rmax = remaind(ifreq)
                        nmax = ifreq
                     endif
 30               continue
c     Add one more frequency corresponding to that mode with the
c     largest remainder (see RN II, 3/30/00)
                  VIB(jcount) = freq(iWell, nmax)
                  ntotal = ntotal + 1
               endif
            endif

c***********************************************************************
c     If frequencies and rotors have been read from file, set this up;
c     because other parts of code use degen vector, this must be set ...
c***********************************************************************
            if(FrRtRd(iWell)) then
               do 32 jcount = 1, nFreqs(iWell)
                  VIB(jcount) = freq(iWell,jcount)
                  degen(iWell,jcount) = 1.0
 32            continue
               sumdegen = dfloat(nFreqs(iWell))
               nfrot = n1drot(iWell)
               ntotal = nFreqs(iWell)

            endif

c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** START                                     c
c----------------------------------------------------------------------c
c RRKM substitution lines -- use same as for well ...
c skip setting of maxBSindex if we are on the 'TS' -- last well and
c RRKM is set ...
               if(RRKM.and.(iWell.eq.ndos)) then
                  maxBSindex(iWell) = maxBSindex(i1RRKM)
                  goto 31
               endif
c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** END                                       c
c----------------------------------------------------------------------c

c***********************************************************************
c     Calculation of max energy level for SBS
c
c     1.  Send along B = 1.0 for the 1d rotors since these technically
c     don't matter (they should cancel)  1.0*dE is added to agree with
c     end of compchem function in cdchem.f
c
c     2.  For chemical activation:  go to 4.5 kT above entrance channel
c     ,or 1.5kT above eMax when 4.5kT above entrance does not exceed
c     eMax + 1.5kT
c
c     3.  For dissoc:  go to 4.5 kT above well bottom; or 1.5 kT above
c     highest exit eMax if 4.5kT does not exceed eMax + 1.5 kT
c
c     Note:  AYC is Modular Man!  No special handling needed for dissoc.
c     Emax will be from bottom of well, and Ewell(inpWell) = 0 for
c     dissoc case; everything works out the same

            xnatoms = (sumdegen + 6.0D0 + dfloat(nrot))/(3.0)

c     eDist = 2.0*(1.5*xnatoms*Rcm*temp) + eMax

c     dmm 20020429  After some problems running Preeti cases in XMG,
c     changed this to be 1/2 of eDist above eMax for criteria (Preeti
c     calc ends up exceeding nmax or maxBSindex)

            eDist = 6.0*(1.5*xnatoms*Rcm*temp)

            if ( eDist.ge.(eMax + 0.5*eDist)) then
               nmax = nint((eDist - Ewell(iWell)/2.859D-3)/bsgs) + 1
            else
               nmax = nint( (eMax + 0.5*eDist-Ewell(iWell)/2.859D-3)
     $              /bsgs)+ 1
            endif
            maxBSindex(iWell) = nmax

c***********************************************************************
c     BARKER OUTPUT Setting of nmax ...
c***********************************************************************

c***********************************************************************
c     For printing Barker output, need nmax to be large enough to agree
c     with 'barker.inp' file ...
            if(barker) then
               open(111, file='barker.inp', status='old')
               read(111,*) Egrain1, imax1, Emax2, Isize, Viblo
               close(111)

c     Notes on 'nmax' for Barker code generation:  Barker code
c     multiwell-1.01 computes k(E) with E = 0 defined as TS ZPE.  This
c     means if calculation goes up to Emax2, then max rho necessary will
c     actually be Emax2 + Ebarmax (maximum barrier out of this well).
c     It's dumb since the actual wells are only calculated to Emax2, but
c     the Barker code (foolishly) follows this format with E=0 at TS
c     ZPE for the k(E) values as well.  These super-high k(E)'s can't
c     actually be used by multiwell (because the rho's only go up to
c     Emax, not Emax +Ebarmax!) but MW barfs w/o them.
c
c     Decided not to rewrite Barker code to fix this inefficiency, so
c     here lies much wasted CPU time until Barker fixes it himself ...
c
c     Note small chance that nmax calc'd this way will not be large
c     enough for good CHEMACT or DISSOC calc, but in practice probably
c     OK.

c***********************************************************************
c     Loops to find max barrier in kcal/mol
c***********************************************************************

               Ebarmax = 0.0
c     First check product channels ...
               do 300 ip1 = 1, nProds(iWell)
                  if(EProd(iWell, ip1).gt.Ebarmax) then
                     Ebarmax = EProd(iWell, ip1)
                  endif
 300           continue

c     Now check isomer channels ...
               do 301 iw1 = 1, mxWells
                  if(EIsom(iWell, iw1).gt.Ebarmax) then
                     Ebarmax = Eisom(iWell, iw1)
                  endif
 301           continue

c***********************************************************************
c     End loops to find max barrier ... convert Ebarmax to kcals!
c***********************************************************************

               Emax2 = Emax2 + Ebarmax*349.775
               nmax = nint(Emax2/bsgs) + nint(2.0*(Emax2/dfloat(Isize))
     $              /bsgs) + 100

               maxBSindex(iWell) = nmax

            endif

c***********************************************************************
c     End BARKER OUTPUT setting of nmax
c***********************************************************************

c***********************************************************************

c     For printing code: find the max nmax
            if (nmax.gt.maxnmax) then
               maxnmax = nmax
            endif
c     Error checking
            if(nmax.gt.maxNOS) then
               write(*,*)
     $   ' ERROR::genRhos.f: Number of grains for sum of states greater'
               write(*,*)
     $   '                   than maxNOS.  Choose larger grainsize or '
               write(*,*)
     $   '                   increase maxNOS in straightBS.fh'

               write(*,*) 'maxNOS = ', maxNOS, ' nmax = ', nmax
               write(*,*) ' ... exiting ...'
               stop
            endif

c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** START                                     c
c----------------------------------------------------------------------c
c     re-entry point below for case where iWell is the RRKM TS

 31         continue
c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** END                                       c
c----------------------------------------------------------------------c

c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** START                                     c
c----------------------------------------------------------------------c
c     Set the rotational constants BR(1) for the k-modes if these have
c     been specified (i.e. krot = 1)
            if(RRKM.and.(rotk(iWell).gt.0.0)) then
               BR(1) = rotk(iWell)
               BDG(1) = 1.0
            else
               BR(1) = 20.0
               BDG(1) = 0.0
            endif
c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** END                                       c
c----------------------------------------------------------------------c

            call SBS(nmax, ntotal, VIB, (nRot + nfrot), BR,BDG,bsgs)
c     Copy results to xrhoofn
            do 60 i=1, nmax
               xrhoofn(i, iWell) = TBS(i)
 60         continue

         endif
c***********************************************************************
c     End if-statement block for recalc of well, continue on to next
c     well
c***********************************************************************

 10   continue

c***********************************************************************
c Printing results for testing (comment out for normal usage)
c
c$$$      open(100, FILE='genRho.out')
c$$$      do 110 i=1, maxnmax
c$$$         write(100, 120) dfloat(i-1)*bsgs, (xrhoofn(i,iWell), iWell=1,
c$$$     $        nWells)
c$$$
c$$$
c$$$
c$$$
c$$$ 110  continue
c$$$ 120  format(F9.1, 6(1X,5E16.9,/))
c***********************************************************************

c***********************************************************************
c Print results for Barker MEQ density file
c
      if(barker) then
         call PrintBarkerDens(maxnmax)
      endif
c***********************************************************************

      return
      end

c***********************************************************************
c dmm 20000404
c
c     FUNCTION xNOS
c
c     This function evaluates the number of states in a width dE given
c     that xrhoofn arrays have been previously evaluated.
c
c     Note that this algorithm EXPECTS bsgs to be considerably lower than
c     dE.  If it is not, it issues a warning but ploughs ahead, at the
c     user's peril.  The results become a bit unpredictable for bsgs
c     larger than the dE grain, jsut as for master equation cases.
c
c     E   = energy in kcal/mol above ZPE of molecule
c     dE  = width, in kcal/mol
c     iWell = well number
c***********************************************************************

      function xNOS(energy, width, iWell)

      implicit real*8 (a-h, o-z)
      include 'cdparams.fh'
      include 'straightBS.fh'
      include 'cdnos.fh'

      integer nE1, nE2

c     Set xNOS to 0.0 if called with a negative energy.  This is natural
c     when coming from getKofEmatrix (signifies no reaction poss. to
c     that well/product at energy)

      if(energy.lt.0.D0) then
         xNOS = 0.0D0
         return
      endif

c     Convert energy and width into cm-1; convert into grainsize index

c     dmm 20011102  Added offset of 1 to nE1 and nE2 ... probably will
c     not affect most results but this is more accurate as xrhoofn
c     begins as xrhoofn(1) at energy level 0.0 cm-1

      xen = energy*349.755
      xwd = width*349.755
      nE1 = nint(xen/bsgs) + 1
      nE2 = nint( (xen + xwd)/bsgs) + 1
c pey 15aug04
c      write(*,*) 'xnos = ', ne1, ne2, maxBSindex(iWell)
c pey
c     Warning check: issue if dE < 2*bsgs and not calculating Q.
      if((xwd.lt.2.0*bsgs).and.(.not.pfflag)) then
         write(*,10)
         write(*,20)
         write(*,30)
      endif

 10   format(/,'WARNING::genRhos.f:xNOS: width for calculation order')
 20   format('           of bsgs. Reset bsgs to lower value or use ')
 30   format('           larger dE.')
c     Warning if nE2 exceeds max. energy of xrhoofn for this well
      if(nE2.gt.maxBSindex(iWell)) then
         write(*,40)
         write(*,50)
         write(*,60)
         write(*,70), iWell, nE2, maxBSindex(iWell)
c         pause
      endif
 40   format(/,'WARNING::genRhos.f:xNOS: xrhoofn index higher than max')
 50   format('           calculated index for this well.  Results will')
 60   format('           be unpredictable.  Use larger nmax.')
 70   format(/,'iWell = ',I4, ' nE2 = ', I10, ' maxBSindex = ', I10)

c     Return xNOS
      xNOS = xrhoofn(nE2, iWell) - xrhoofn(nE1, iWell)

      return
      end
c dmm 2000404
c***********************************************************************

c***********************************************************************
c***********************************************************************
c     FUNCTION xNOScm
c
c     Same as xNOS but assumes energy in cm-1
c***********************************************************************
c***********************************************************************

      function xNOScm(energy, width, iWell)

      implicit real*8 (a-h, o-z)
      include 'cdparams.fh'
      include 'straightBS.fh'
      include 'cdnos.fh'

      integer nE1, nE2

c     Set xNOS to 0.0 if called with a negative energy.  This is natural
c     when coming from getKofEmatrix (signifies no reaction poss. to
c     that well/product at energy)

      if(energy.lt.0.D0) then
         xNOS = 0.0D0
         return
      endif

c     Convert energy and width into cm-1; convert into grainsize index

c     dmm 20011102  See above for explanation of offset on nE1 and nE2
      xen = energy
      xwd = width
      nE1 = nint(xen/bsgs) + 1
      nE2 = nint( (xen + xwd)/bsgs) + 1

c     Warning check: issue if dE < 2*bsgs and not calculating Q.
      if((xwd.lt.2.0*bsgs).and.(.not.pfflag)) then
         write(*,10)
         write(*,20)
         write(*,30)
      endif

 10   format(/,'WARNING::genRhos.f:xNOS: width for calculation order')
 20   format('           of bsgs. Reset bsgs to lower value or use ')
 30   format('           larger dE.')
c     Warning if nE2 exceeds max. energy of xrhoofn for this well
      if(nE2.gt.maxBSindex(iWell)) then
         write(*,40)
         write(*,50)
         write(*,60)
         write(*,70), iWell, nE2, maxBSindex(iWell)
      endif
 40   format(/,'WARNING::genRhos.f:xNOS: xrhoofn index higher than max')
 50   format('           calculated index for this well.  Results will')
 60   format('           be unpredictable.  Use larger nmax.')
 70   format(/,'iWell = ',I4, ' nE2 = ', I10, ' maxBSindex = ', I10)

c     Return xNOS
      xNOScm = xrhoofn(nE2, iWell) - xrhoofn(nE1, iWell)

      return
      end
c dmm 2000404
c***********************************************************************

c***********************************************************************
c***********************************************************************
c dmm 20000627
c
c     SUBROUTINE InitGenRhos
c
c     All this subroutine does is zero out the wellscalcd vector (see
c     genRhos, cdnos.fh).  It is probably not necessary, but it helps
c     make sure.  Called from chemdis2.f ... it's clunky but it
c     works for now.  A more robust and elegant solution escapes me.
c
c***********************************************************************
c***********************************************************************

      subroutine InitGenRhos

      implicit none
      include 'cdparams.fh'
      include 'straightBS.fh'
      include 'cdnos.fh'

      integer iWell

      do 10 iWell = 1, mxWells
         wellscalcd(iWell) = 0.0D0
 10   continue

      return
      end

c***********************************************************************
c***********************************************************************
c dmm 20011216
c
c     SUBROUTINE WRSumDens
c
c     This function computes the Whitten Rabinovitch approx to the sums
c     and densities of states; note that Albert & orig. authors have
c     their own version.  I don't trust its units consistency with
c     direct count, hence this function.  Plan is for this to be called
c     above the "break" energy of the double-array density of states
c     files produced by PrintBarkerDens for the multiwell code.
c
c     PLACEHOLDER:  possible serious issue at higher T and P by using WR
c     approx for DOS/SOS but 3-freq direct count for QRRK k(E)s.
c     Probably better to use the smoother WR where possible anyway,
c     since 3-freq approx is only meaningful in a "smooth" sense, and
c     spikiness/noise resulting from direct count is artificial.
c
c     Evr:  energy in cm-1 above ZPE
c     dens:  density of states returned in 1/cm-1 (molecular)
c
c***********************************************************************
c***********************************************************************

      SUBROUTINE WRSumDens(iWell, Evr, Sum, Dens)

      implicit none

      integer iWell, i, j, ndg
      real*8 Evr, Sum, Dens, piv, Eprime, Ezero
      real*8 vmsq, vsqm, vfsqsum, vfsum, betaR, xnrot, xnvib, Qrot, Bk,
     $     Bdegen, wfac, afrac, lprodf, lWRsum, gamarg, dwdEp, lWRdens

      double precision gammln

      include 'cdparams.fh'
      include 'cdisprop.fh'

      piv = 3.14159

c***********************************************************************
c     Begin constants evaluation -- these do not change given same
c     molecule
c***********************************************************************

c     Calculate zero-point energy

      Ezero = 0.0
      vfsum = 0.0
      vfsqsum = 0.0
      lprodf = 0.0
      xnvib = 0.0
      do 10 i = 1, nFreqs(iWell)
         ndg = nint(degen(iWell,i))
         do 15 j = 1, ndg
            xnvib = xnvib + 1.0
            Ezero = Ezero + 0.5*freq(iWell,i)
            vfsum = vfsum + freq(iWell,i)
            vfsqsum = vfsqsum + (freq(iWell,i))**(2.0)
            lprodf = lprodf + LOG(freq(iWell,i))
 15      continue
 10   continue

c     Total (float) number of rotational modes
      xnrot = dfloat(nRot + n1drot(iWell))

c     Calculate vmsq, vsqm
      vsqm = vfsqsum/xnvib
      vmsq = (vfsum/xnvib)**(2.0)

c     betaR for the WR approx
      betaR = ((xnvib-1.0)/(xnvib))*((xnvib+0.5*xnrot)/(xnvib))*(
     $    vsqm/vmsq)

c     Evaluate rotational partition function: k-mode
c     Note: rotational constant default here of 20.0 MUST agree with that
c     set above in genRhos!
c     PLACEHOLDER:  there is currently no handling of rotational mode
c     degeneracies.  These are all assumed 1 by the code (should be
c     changed at some point) since QRRK uses only ratios of DOS anyway,
c     where the degeneracy would cancel out.

      Qrot = 1.0
      if(nRot.gt.0) then
         if(rotk(iWell).le.0.0) then
            Bk = 20.0
            Bdegen = 1.0
         else
            Bk = rotk(iWell)
            Bdegen = 1.0
         endif
         Qrot = Qrot*((piv/(Bk*Bdegen))**(0.5))
      endif

c     Evaluate other rotors.  PLACEHOLDER:  Again, rotational modes are never
c     actually read in for internal rotors, nor stored, since the QRRK
c     k(E) formula causes them to cancel along with degeneracies.  This
c     is probably a deficiency which should be corrected.  Constant
c     here should be same as the default set around line 70 of genRhos.f

      if(n1drot(iWell).gt.0) then
         do 20 i = 1, n1drot(iWell)
            Bk = 20.0
            Bdegen = 1.0
            Qrot = Qrot*((piv/(Bk*Bdegen))**(0.5))
 20      continue
      endif

c***********************************************************************
c End constants evaluation
c***********************************************************************

c     Evaluate Eprime
      Eprime = Evr/Ezero
c     Test function limits
      if(Eprime.lt.0.1) then
         write(*,900)
         stop
      else if(Eprime.gt.8.0) then
         write(*,901)
      endif
c     Calculate wfac
      if((Eprime.lt.1.0)) then
         wfac = (5.0*Eprime + 2.73*(Eprime**0.5) + 3.51)**(-1.0)
         dwdEp = -(5.0 + 1.365*(Eprime**(-0.5)))*(wfac**2.0)
      else
         wfac = EXP(-2.4191*(Eprime**0.25))
         dwdEp = -(0.60478*(Eprime**(-0.75)))*wfac
      endif
c     Calculate afrac
      afrac = 1.0 - betaR*wfac
c     Calculate gammln function argument
      gamarg = xnvib + 0.5*xnrot + 1.0
c     Calculate log of WR sum of states and dens of states
      lWRsum = LOG(Qrot) - gammln(gamarg) + (gamarg-1.0)*LOG(Evr + afrac
     $     *Ezero) - lprodf
      Sum = EXP(lWRsum)

      gamarg = xnvib+0.5*xnrot
      lWRdens = (LOG(Qrot) - gammln(gamarg) + (gamarg-1.0)*LOG(Evr +
     $     afrac*Ezero) - lprodf) + LOG(1.0 - betaR*dwdEp)

      Dens = EXP(lWRdens)

 900  format(1x,'ERROR:  WR called with energy too low for accuracy',1x,
     $    ' ... stopping')
 901  format(1x
     $     ,'WARNING:  WR called with energy greater than 8 times ZPE')

      return
      end

c***********************************************************************
c***********************************************************************

c $Id$
c $Author$
c $Date$
c $Source$
c $Revision$
c $Log$
c Revision 1.2  2007-02-21 03:16:25  sandeeps
c Changed a minor syntax so that g95 can compile it. No major change in the program itself
c
c Revision 1.1.1.1  2007/02/20 23:10:23  sandeeps
c Fortran_software contain
c 1. therfir
c 2. chemdis
c 3. reactorModel (chemkin, as modified by Paul)
c 4. ODESolver (the fortran code which RMG uses to call daspk or dassl)
c 5. fit3pbnd
c
c Revision 1.16  2001/12/16 20:44:45  dmmatheu
c before adding Whitten-Rabinovitch subroutine, to get smooth DOS/SOS
c profile above high energy for Barker file input ...
c
c Revision 1.15  2001/12/02 17:51:06  dmmatheu
c Before offset fix:  it seems that nE1 and nE2 should be incremented by
c 1 automatically, since xrhoofn(1) is the same as energy level 0.  Not
c sure how this error got past for so long.  This gives "look-ahead" dos
c in xNOScm:  call at energy, width gives number of states from energy
c to energy+width.
c
c Revision 1.14  2001/10/19 16:02:04  dmmatheu
c before changes to generate SOS for the 1 RRKM rate constant
c
c Revision 1.13  2001/06/13 14:12:16  dmmatheu
c After changes to turn Barker file outputs on/off, major debug of
c mxWells and mxProds (zeroing loops of MakeShadow were exceeding array
c bounds).  Also added xNOScm from ASA-MasterEqn.. directory ... merging
c to one chemdis file
c
c Revision 1.12  2000/08/09 18:37:31  dmmatheu
c after tlada bug fix on dissoc
c
c Revision 1.11  2000/06/27 22:12:16  dmmatheu
c Changed to avoid recalculating density of states for a given well,
c when it's not necessary.  Corresponding changes to chemdis2.f,
c cdnos.fh.  Checked on standard vinyl+O2 9P14T example, looks good (no
c ASA attempted).
c
c Revision 1.10  2000/06/23 17:54:34  dmmatheu
c before Tlada merge ...
c
c Revision 1.9  2000/06/20 18:47:28  dmmatheu
c Bug in previous version never discovered ... things seem to be working
c fine now (??).  Very strange.
c
c Revision 1.8  2000/06/20 18:18:47  dmmatheu
c Before working on very strange bug -- commenting out the genRho.out
c file production lines causes different results in fort.65 for dissoc
c case
c
c Revision 1.7  2000/05/09 18:05:30  dmmatheu
c Works for vinyl + O2 example from May 4.
c Before changes with dissoc implementation.
c
c Revision 1.6  2000/04/04 21:13:50  dmmatheu
c Applied 3kT rule to estimate max for density of states ... need to
c build safety AND warning for exceeding this into code ...
c
c Revision 1.5  2000/04/04 20:38:51  dmmatheu
c Code compiles and runs without core-dumps.  eMax/max height for genRho
c is not handled correctly.  Need to estimate eMax.
c
c Revision 1.4  2000/04/04 15:12:04  dmmatheu
c Check in before adding xNOS function
c
c Revision 1.3  2000/03/30 23:43:18  dmmatheu
c Successful implementation in CHEMDIS on vinyl + O2 (generation of
c xrhoofn only)
c
c Revision 1.2  2000/03/30 20:32:23  dmmatheu
c pre-compiled version
c

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


! Copyright (c) 2012-2013 Bartosz Szczesny
! LICENSE: The MIT License (MIT)

program equivalent_variables
    integer :: i, j, n( 3 ), m( 3 )
    equivalence ( i, j )
    equivalence ( n, m )

    i = 1
    print "( a, i1, a, i1 )", "i = ", i, ", j = ", j

    j = 2
    print "( a, i1, a, i1 )", "i = ", i, ", j = ", j

    n( 1 ) = 11
    n( 2 ) = 22
    n( 3 ) = 33
    print "( a, 3i3 )", "m =", m

    stop
end


! @Copyright 2007 Kristjan Haule

SUBROUTINE Cmp_Optics(fenergy, fUdmft, fsymop, fh_m, fhb, ommax, Temperature_, delta, Nd, gammac, gamma, alphaV, Ndirection, Qsym, Qsimple, nom, nat, iso, norbitals, ncix, natom, nkpt, nmat, nume, Qcomplex, lmax2, maxdim2, maxdim, maxsize, dwindow, InterbandOnly)
  USE com_mpi, ONLY: nprocs, myrank, master, Reduce_MPI, Reduce_MPI_dos
  IMPLICIT NONE
  CHARACTER*100, intent(in) :: fenergy, fUdmft, fsymop!, fmommat
  REAL*8, intent(in)  :: ommax, delta, gammac, gamma, alphaV(3,3,Ndirection), dwindow, Temperature_
  INTEGER, intent(in) :: Nd, Ndirection, fh_m
  LOGICAL, intent(in) :: Qsym, Qsimple
  INTEGER, intent(in) :: fhb, nom, nat, iso, norbitals, ncix, natom, nkpt, nmat, nume
  LOGICAL, intent(in) :: Qcomplex
  INTEGER, intent(in) :: lmax2, maxdim2, maxdim, maxsize
  LOGICAL, intent(in) :: InterbandOnly
  ! Functions
  !COMPLEX*16 :: Opt_Inside_Simple, Opt_Inside, Opt_Noncorrelated
  interface
     REAL*8 Function ferm(x)
       REAL*8, intent(in):: x
     end function ferm
  end interface
  ! locals
  COMPLEX*16 :: ctmp(Ndirection), ctmp2(Ndirection), copt(Ndirection)
  CHARACTER*10:: skii
  LOGICAL :: Tcompute, pform
  INTEGER :: ios, pr_proc
  INTEGER :: itape, i, j, k, l, is, N, nemin, nemax, ikp, iikp, ip, iq, NE, nb_min, nb_max, dir, n_min, n_max
  REAL*8  :: EF, VOL, EMIST, renorm_wgh, wommax, wdelta, wgamma, wdwindow, Ek(nume), wgh(nkpt), womega, eps_p, eps_m, ferm_factors, Temperature
  integer :: fh_p, fh_s, fh_d, fh_o!, fh_m
  integer :: numk, nsymop, nbands, isym, ddi, iw, istart, iend, ix, iw_p, iw_m
  integer :: iorb, wkii, wkis, iom, iband, nord, N0, Nw, Nsymw, icix, tnorbitals
  integer :: nindo(norbitals), csize(ncix), tnindo(norbitals), cixdim(ncix), nl(natom), ll(natom,4)
  INTEGER :: cix(natom,4), iorbital(natom,4), nind(natom,4), iSx(maxdim2,norbitals), Sigind(maxdim,maxdim,ncix)
  INTEGER,    ALLOCATABLE :: amesh(:), bandind(:)
  REAL*8,     ALLOCATABLE :: zomega(:), conduc(:,:), opimat(:,:,:), omega(:)
  COMPLEX*16, ALLOCATABLE :: STrans(:,:,:,:), DMFTU(:,:,:), gc(:)
  COMPLEX*16, ALLOCATABLE :: sigma(:,:,:), wsigma(:,:), gij(:,:), zekw(:,:), Al(:,:,:), Ar(:,:,:)
  COMPLEX*16, ALLOCATABLE :: Cv(:,:,:), Dv(:,:,:), Vel(:,:,:), Ve(:,:,:), Cv0(:,:,:)
  COMPLEX*16, ALLOCATABLE :: L1(:,:,:), L2(:,:,:), L3(:,:,:), L4(:,:,:), Ltmp(:,:)
  COMPLEX*16, ALLOCATABLE :: lg_q(:), lg_p(:)
  LOGICAL,    ALLOCATABLE :: Qinside(:)
  ! constants
  REAL*8, PARAMETER      :: Ry2eV = 13.60569193
  REAL*8, PARAMETER      :: elresunit = 2.173985e-5  ! [ aBohr*hbar/ecoul^2 ] in units of Ohm*cm !!! Natural unit for electrical resistivity: [ aBohr*hbar/ecoul^2 ]
  REAL*8, PARAMETER      :: PI = 3.14159265358979
  COMPLEX*16, PARAMETER  :: IMAG = (0.0D0,1.0D0)
  REAL*8 :: COmega(2*Nd+1)

  ! Converting to Rydbergs
  wommax = ommax/Ry2eV
  wdelta = delta/Ry2eV
  wgamma = gamma/Ry2eV
  wdwindow = dwindow/Ry2eV
  Temperature = Temperature_/Ry2eV

  fh_p = 399 ! file for DMFT transformation UDMFT

  ! Logarithmic frequency mesh for integration
  N0 = int(log((wommax/wdelta)/Nd)/log(2.) + 1 + 0.5) ! Number of different logarithmic intervals
  if (N0 .LT. 1) N0=1 ! should not be zero
  Nw = (N0+1)*Nd/2   ! Number of all points in the logarithmic mesh
  ALLOCATE ( amesh(2*Nd+1), zomega(2*Nw), conduc(Nw,Ndirection), gc(2*Nw) )
  gc=0
  print *, 'Number of logarithmic intervals=', N0, 'Number of all points in the logarithmic mesh=', Nw

  ! Many index arrays are imported from DMFT
  CALL Read_The_Rest_Basic_Arrays(fhb, nindo, cixdim, nl, ll, iorbital, cix, nind, csize, iSx, Sigind, EF, VOL, norbitals, ncix, natom, maxdim, maxdim2)

  ! Reading the self-energy for all orbitals
  ALLOCATE( sigma(nom,maxsize,ncix), omega(nom), wsigma(maxsize,ncix) )
  do icix=1,ncix
     itape = 80+icix
     WRITE(skii,fmt='(I2)') icix
     open(itape, file=TRIM('sig.inp')//ADJUSTL(TRIM(skii)), status='old', form='formatted')
     CALL ReadSelfenergy(itape, sigma(:,:,icix), omega, gammac, csize(icix), nom, maxsize)
  enddo
  ! Self-energy was in eV -> transform to Ry
  omega(:) = omega(:)/Ry2eV
  sigma(:,:,:) = sigma(:,:,:)/Ry2eV

  ! Kohn-Sham eigenvalues from the file
  open (59, file=fenergy, status='old', form='formatted')

  ! Reading symmetrization for matrix elements. It is used to transform from an irreducible k-point to general k-point.
  fh_s = 98
  open(fh_s, file=fsymop, status='old')
  READ(fh_s, fmt='(I6)') nord
  ALLOCATE( opimat(3,3,nord) )
  DO k=1,nord
     READ(fh_s, fmt='(3(3f8.5/))') ((opimat(j,i,k),i=1,3),j=1,3)
  END DO

  !! Reading linearization energies in Energy file
  itape=59
  do i=1,nat
     READ(itape, fmt='(f9.5)') EMIST !---- At the beginninge we have linearization energies --!
     READ(itape, fmt='(f9.5)') EMIST !---- We just skip them ---------------------------------!
  enddo

  wkis=0 ! DMFT transofmration might be in many files. Starting by the first.
  ! DMFT transformation
  WRITE(skii,fmt='(I2)') wkis

  !pform = .true.
  pform = .false.
  if (pform) then
     open(fh_p, file=TRIM(fUdmft)//ADJUSTL(TRIM(skii)), status='old',form='formatted')
  else
     open(fh_p, file=TRIM(fUdmft)//ADJUSTL(TRIM(skii)), status='old',form='unformatted')
  endif

  if (pform) then
     READ(fh_p,*) numk, nsymop, tnorbitals
     READ(fh_p,*) (tnindo(iorb), iorb=1,norbitals)
  else
     READ(fh_p) numk, nsymop, tnorbitals
     READ(fh_p) (tnindo(iorb), iorb=1,norbitals)
  endif

  ! optics matrix elements
  !open(fh_m, file=fmommat, status='old')
  READ(fh_m,*)

  COmega=0

!!! ---------- Preparation of arrays for paralel executaion --------------
  pr_proc  = floor(nkpt/DBLE(nprocs)+0.999)  ! The maximum number of points calculated per processor
  WRITE(6,'(A,I4,2x,A,I4)') 'pr_proc=', pr_proc, 'tot-k=', nkpt

  wgh=0
  conduc = 0
  wkii=0
  DO ikp=1,nkpt   ! over all irreducible k-points
     !--- We need to go over all k-points even though we will compute only some of them on this processor.
     !--- This is because we need to read vector file sequentially.
     iikp = ikp-myrank*pr_proc            ! The index of the point to compute. If negative, do not compute!
     if (iikp.GT.pr_proc) EXIT            ! Processor finished. The rest of the points will be taken care of by the other processors.
     Tcompute=.FALSE.
     if (iikp.gt.0) Tcompute=.TRUE.       ! If Tcompute is true, the point needs to be computed.
     ! Reading band energies from energy file
     CALL ReadEnergiesK(Ek, wgh, NE, ikp, nkpt, nume)
     ! Start reading DMFT transformation
     CALL ReadDMFT_TransK_Outside(nsymop, nbands, nemin, numk, wkii, wkis, fUdmft, fh_p, pform, ikp)
     nemax = nemin+nbands-1

     ALLOCATE( bandind(NE) )

     DO ip=1,NE
        bandind(ip)=ip
     ENDDO
     DO ip=1,NE
        do iq=ip+1,NE
           if (abs(Ek(ip)-Ek(iq)).LT.1e-11) then
              bandind(iq)=bandind(ip)
           endif
        enddo
     ENDDO

     ALLOCATE( Qinside(NE) )
     Qinside=.FALSE.
     do ip=1,nbands
        Qinside(ip+nemin-1) = .TRUE.
     enddo

     ALLOCATE ( DMFTU(nbands,maxdim2,norbitals) )
     ALLOCATE ( Vel(NE,NE,3) )
     IF (Tcompute) THEN
        ALLOCATE( STrans(maxsize,ncix,nbands,nbands) )
        allocate( gij(nbands,nbands), zekw(NE,2*Nd+1), Al(nbands,nbands,2*Nd+1), Ar(nbands,nbands,2*Nd+1))
        ALLOCATE( Ve(NE,NE,3), Cv0(NE,NE,Ndirection) )
        ALLOCATE( Cv(NE,NE,Ndirection), Dv(NE,NE,Ndirection) )
        ALLOCATE( L1(NE,NE,3), L2(NE,NE,3), L3(NE,NE,3), L4(NE,NE,3), Ltmp(nbands,nbands) )
        ALLOCATE( lg_q(NE), lg_p(NE) )
     ELSE
        wgh(ikp)=0
     ENDIF

     CALL ReadVelocityK(Vel, NE, fh_m, nb_min, nb_max) !nemin, nbands)

     if (Qsym) then
        Nsymw = nsymop  ! Over all reducible k-points (all group transformations)
     else
        Nsymw = 1       ! Only irreducible k-points needed
     endif

     DO isym=1, nsymop

        ! Continuing reading DMFT transformation
        CALL ReadDMFT_TransK_Inside(DMFTU, fh_p, pform, isym, nindo, nbands, norbitals, maxdim2)

        if (isym.GT.Nsymw .OR. .NOT.Tcompute) CYCLE ! If only the average optics is needed, we can go over the irreducible BZ only.

        !print *, ikp, isym

        ! For more efficient transformation of self-energy, we create special array of transformation STrans
        CALL CompressSigmaTransformation2(STrans, DMFTU, Sigind, iSx, cix, csize, iorbital, ll, nl, natom, iso, ncix, maxdim, maxdim2, lmax2, norbitals, nbands, maxsize)

        ! Changing matrix elements to general k-point
        CALL Transform_To_Reducible(Ve, isym, Vel, opimat, NE, nord)

        !!!! For non-correlated
        L1 = Ve
        L2 = Ve
        L3 = Ve
        L4 = Ve
        DO ip=nb_min,nb_max
           DO iq=nb_min,nb_max
              ctmp=0
              DO dir=1,Ndirection
                 DO j=1,3
                    DO l=1,3
                       ctmp(dir) = ctmp(dir) + Ve(ip,iq,j)*Ve(iq,ip,l)*alphaV(l,j,dir)
                    ENDDO
                 ENDDO
              ENDDO
              Cv0(ip,iq,:) = ctmp
              Cv(ip,iq,:) = ctmp
              Dv(ip,iq,:) = ctmp
           ENDDO
        ENDDO
        ! For non-correlated
        do i=1,2*Nd+1
           zekw(:,i) = Ek(:)
        enddo

        iw=1
        ddi = 1
        do j=1,N0
           ! The algoritm for the frequency mesh is the following:
           ! We start with Nd negative and Nd positive points of equidistant mesh. Incluing zero, we have 2*Nd+1 points.
           ! The points are [-Nd,-Nd+1,...0, 1,.., Nd]*wdelta
           ! With this mesh, we can compute all integrals up to frequency Nd*wdelta, where wdelta is the distance between the points
           ! In the second step, we keep every second point: [-Nd,-Nd+2,...0,2,...Nd]*wdelta
           ! and we add additional points, so that we have again 2*Nd+1 points. We thus have the following points
           ! [-2*Nd, -2*Nd+2,...0,2,..2*Nd]*wdelta.
           ! The old points (one half of all the points) do not need to be recomputed. Only the new points are recomputed.
           ! The interval is doubles N0 times. N0 depends on the cutoff ommax.
           if (j.EQ.1) then
              ix=1
              do i=-Nd,Nd
                 amesh(i+Nd+1) = ddi*i
                 CALL FindEigsys(zekw(nemin:nemax,i+Nd+1), Al(:,:,i+Nd+1), Ar(:,:,i+Nd+1), ix, wsigma, gij, amesh(i+Nd+1)*wdelta, omega, STrans, Ek, nom, ncix, nbands, nemin, maxsize, nume, csize, sigma)
              enddo
           else
              do i=1,Nd/2 ! These points were computed in previous step. Just reshuffle them into the right place.
                 amesh( i+Nd+1) = amesh( 2*i+Nd+1)
                 zekw(nemin:nemax,i+Nd+1) = zekw(nemin:nemax,2*i+Nd+1)
                 Al(:,:,i+Nd+1) = Al(:,:,2*i+Nd+1)
                 Ar(:,:,i+Nd+1) = Ar(:,:,2*i+Nd+1)
              enddo
              do i=-1,-Nd/2,-1 ! These points were computed in previous step. Just reshuffle them into the right place.
                 amesh( i+Nd+1) = amesh( 2*i+Nd+1)
                 zekw(nemin:nemax,i+Nd+1) = zekw(nemin:nemax,2*i+Nd+1)
                 Al(:,:,i+Nd+1) = Al(:,:,2*i+Nd+1)
                 Ar(:,:,i+Nd+1) = Ar(:,:,2*i+Nd+1)
              enddo
              do i=-Nd,-Nd/2-1 ! These points are new. Need to compute eigensystem for them.
                 amesh(i+Nd+1) = ddi*i
                 CALL FindEigsys(zekw(nemin:nemax,i+Nd+1), Al(:,:,i+Nd+1), Ar(:,:,i+Nd+1), ix, wsigma, gij, amesh(i+Nd+1)*wdelta, omega, STrans, Ek, nom, ncix, nbands, nemin, maxsize, nume, csize, sigma)
              enddo
              do i=Nd/2+1,Nd,1 ! These points are new. Need to compute eigensystem for them.
                 amesh(i+Nd+1) = ddi*i
                 CALL FindEigsys(zekw(nemin:nemax,i+Nd+1), Al(:,:,i+Nd+1), Ar(:,:,i+Nd+1), ix, wsigma, gij, amesh(i+Nd+1)*wdelta, omega, STrans, Ek, nom, ncix, nbands, nemin, maxsize, nume, csize, sigma)
              enddo
           endif

           if (j.EQ.1) then
              istart=1
              iend=Nd
           else
              istart=Nd/2+1
              iend=Nd
           endif
           if (Temperature.NE.0.0) then
              istart=1
              iend=Nd-4*Temperature/wdelta
              if (iend<1) iend=1
           endif
           WRITE(*,'(A,I4,1x,A,I3,1x,A,I3,1x,A,I3,1x,A,I4,1x,A,I4,1x,A,F10.3)') 'ikp=', ikp, 'isym=', isym, 'j=', j, 'istart=', istart, 'iend=', iend, 'Nd=', Nd, '4*T/d=', 4*Temperature/wdelta

           ix=1
           do k=istart,iend
              ! Frequency mesh for printing
              zomega(Nw+iw)   = amesh(k+Nd+1)*wdelta ! positive frequencies
              zomega(Nw-iw+1) = amesh(Nd-k+1)*wdelta ! negative frequencies
              ! For Density of state
              do ip=1,nbands
                 iband = ip+nemin-1
                 gc(Nw+iw)   = gc(Nw+iw)   + wgh(ikp)/(zomega(Nw+iw)  +EF-zekw(iband,Nd+k+1)+wgamma*IMAG)
                 gc(Nw-iw+1) = gc(Nw-iw+1) + wgh(ikp)/(zomega(Nw-iw+1)+EF-zekw(iband,Nd-k+1)+wgamma*IMAG)
              enddo

              womega = amesh(k+Nd+1)*wdelta

              n_min = nb_min
              n_max = nb_max
              do ip=nb_min,nb_max
                 if ( (Ek(ip)-EF).LT.(-womega-wdwindow) ) n_min = ip
              enddo
              do ip=nb_max,nb_min,-1
                 if ( (Ek(ip)-EF).GT.(womega+wdwindow) ) n_max = ip
              enddo

              ctmp=0  ! integration over epsilon=[0,omega]
              if (Temperature.NE.0.0) then !! FINITE TEMPERATURE
                 ! !$OMP PARALLEL DO SHARED(COmega,amesh,Al,Ar,Ve,alphaV) PRIVATE(iw_p,iw_m,eps_p,eps_m,ferm_factors,Cv,Dv,L1,L2,L3,L4,Ltmp,copt) SCHEDULE(STATIC) &
                 ! !$OMP& REDUCTION(+:ctmp)
                 do l=k+1-Nd,Nd+1
                    iw_p = Nd+l                 ! eps
                    iw_m = Nd-k+l               ! eps-omega

                    eps_p = amesh(iw_p)*wdelta    ! eps
                    eps_m = amesh(iw_m)*wdelta    ! eps-omega
                    ferm_factors = ferm(eps_m/Temperature)-ferm(eps_p/Temperature)

                    CALL cmp_C_and_D(Cv, Dv, Al(:,:,iw_p), Ar(:,:,iw_p), Al(:,:,iw_m), Ar(:,:,iw_m), Ve, alphaV, Ndirection, L1, L2, L3, L4, Ltmp, Qinside, nbands, NE, nemin, nemax, n_min, n_max)

                    if (Qsimple) then
                       CALL Opt_Inside_Simple(copt, Cv, Dv, zekw, k, -10000, amesh, iw_p, iw_m, wdelta, wgamma, EF, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, n_min, n_max, InterbandOnly)
                       ctmp = ctmp + copt*ferm_factors
                    else
                       CALL Opt_Inside(copt, Cv, Dv, zekw, k, -10000, amesh, iw_p, iw_m, wdelta, wgamma, EF, lg_q, lg_p, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, n_min, n_max, InterbandOnly)
                       ctmp = ctmp + copt*ferm_factors
                    endif

                    !if (k.eq.1) WRITE(*,'(A,I3,1x,A,F12.6,1x,A,I3,1x,A,F12.6,1x,A,F12.6,1x,A,F12.6)') 'iw_p=', iw_p, 'eps_p=', eps_p/Temperature, 'iw_m=', iw_m, 'eps_m=', eps_m/Temperature, 'ferm=', ferm_factors, 'opt=', -dble(copt(1))
                    if (k.eq.1) then
                       COmega(iw_p) = COmega(iw_p) - dreal(copt(1))*wgh(ikp)*(4*pi)/(2*pi**2*elresunit*VOL)
                    endif
                 enddo
                 ! !$OMP END PARALLEL DO
              else
                 do l=1,k+1
                    iw_p = Nd+l                 ! eps
                    iw_m = Nd-k+l               ! eps-omega

                    CALL cmp_C_and_D(Cv, Dv, Al(:,:,iw_p), Ar(:,:,iw_p), Al(:,:,iw_m), Ar(:,:,iw_m), Ve, alphaV, Ndirection, L1, L2, L3, L4, Ltmp, Qinside, nbands, NE, nemin, nemax, n_min, n_max)

                    if (Qsimple) then
                       CALL Opt_Inside_Simple(copt, Cv, Dv, zekw, k, l, amesh, iw_p, iw_m, wdelta, wgamma, EF, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, n_min, n_max, InterbandOnly)
                       ctmp = ctmp + copt
                    else
                       CALL Opt_Inside(copt, Cv, Dv, zekw, k, l, amesh, iw_p, iw_m, wdelta, wgamma, EF, lg_q, lg_p, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, n_min, n_max, InterbandOnly)
                       ctmp = ctmp + copt
                    endif
                 enddo
              endif

              !!! For non-correlated
              !womega = amesh(k+Nd+1)*wdelta ! omega
              CALL Opt_Noncorrelated(copt,Cv0,Ek,womega,EF,wgamma,nb_min,nb_max,nemin,nemax,nume,Qinside,NE,Ndirection,InterbandOnly)
              ctmp2 = copt

              DO dir=1,Ndirection
                 ! We add proper dimension to optics
                 conduc(iw,dir) = conduc(iw,dir) - dreal(ctmp(dir) + ctmp2(dir))*wgh(ikp)*(4*pi)/(2*pi**2*elresunit*VOL)
              ENDDO
              iw = iw + 1
           enddo
           ddi = ddi * 2
        enddo
     ENDDO
     DEALLOCATE( bandind )
     DEALLOCATE( Qinside )
     DEALLOCATE( DMFTU, Vel )
     IF (Tcompute) THEN
        DEALLOCATE( lg_q, lg_p )
        DEALLOCATE( L1, L2, L3, L4, Ltmp )
        DEALLOCATE( Cv, Dv, Ve, Cv0 )
        DEALLOCATE( gij, zekw, Al, Ar)
        DEALLOCATE( STrans )
     ENDIF
  ENDDO

  ! k-points weights are not normalized. Need to normalize at the end!
  renorm_wgh=0
  DO ikp=1,nkpt
     renorm_wgh = renorm_wgh + wgh(ikp)
  ENDDO

  CALL Reduce_MPI(conduc, COmega, Nw, Ndirection, renorm_wgh, Nd)
  CALL Reduce_MPI_dos(gc,Nw)

  IF (myrank.EQ.master) THEN

     ! Normalization over irreducible k-points and over symmetrization (all k-points)
     gc = gc/(renorm_wgh*Nsymw)
     conduc = conduc/(renorm_wgh*Nsymw) * (2.0/iso)  ! changed in 2019. Needs to be multiplied by 2 due to spin, except if spin-orbit is included

     ! Output of total DOS and optics
     fh_d = 501
     open(fh_d, file='optdos.dat', status='unknown', form='formatted')
     WRITE(fh_d,*) '# omega[eV]   DOS[1/eV]'
     DO iom=1,2*Nw
        if (zomega(iom).NE.0) WRITE(fh_d,*) zomega(iom)*Ry2eV, -aimag(gc(iom))/PI/Ry2eV
     ENDDO
     close(fh_d)

     fh_o = 502
     open(fh_o, file='optics.dat', status='unknown', form='formatted')
     WRITE(fh_o,*) '# omega[eV]   sigma[1/(Ohm*cm)]'
     DO iom=1,Nw
        if (zomega(Nw+iom).GT.0.0) then
           WRITE(fh_o, '(f14.10,2x)', advance='no') zomega(Nw+iom)*Ry2eV
           DO dir=1,Ndirection
              WRITE(fh_o, '(f20.5,1x)', advance='no') conduc(iom,dir)/(zomega(Nw+iom))
           ENDDO
           WRITE(fh_o,*)
        endif
     ENDDO

     if (Temperature.NE.0.0) then
        COmega = COmega/(renorm_wgh*Nsymw)
        fh_d = 503
        open(fh_d, file='what_we_integrate.dat', status='unknown', form='formatted')
        WRITE(fh_d,*) '# omega[eV]   OMEGA[1/eV]'
        DO iom=2,2*Nd+1
           WRITE(fh_d,*) (iom-Nd-1)*wdelta*Ry2eV, COmega(iom)/wdelta
        ENDDO
        close(fh_d)
     endif

  ENDIF

998 CONTINUE

  close(fh_s)
  close(fh_p)
  close(59)
  close(60)
  DEALLOCATE( amesh, zomega, conduc, gc )
  DEALLOCATE( opimat )
  DEALLOCATE( sigma, omega, wsigma )
9010 FORMAT(/,2X,' KP:',I6,' NEMIN NEMAX : ',2I5, ' dE:',2f5.2,' K:',a10 /)
9040 FORMAT(3X,2I4,6E13.6,F13.8)
END SUBROUTINE Cmp_Optics

SUBROUTINE ReadVelocityK(Vel, NE, fh_m, nb_min, nb_max)!, nemin, nbands)
  IMPLICIT NONE
  COMPLEX*16, intent(out) :: Vel(NE,NE,3)
  INTEGER, intent(out) :: nb_min, nb_max
  INTEGER, intent(in) :: fh_m, NE!, nemin, nbands
  ! locals
  INTEGER      :: iikp, i1, i2, ib1, ib2, ixyz
  REAL*8       :: wEmin, wEmax, dE
  CHARACTER*10 :: kname
  COMPLEX*16   :: O(3)

  Vel = 0
!!! Reading optics matrix elements
  READ(fh_m,9010) iikp,nb_min,nb_max,wEmin,wEmax,kname
  do i1=nb_min,nb_max
     do i2=i1,nb_max
        READ(fh_m, 9040) ib1,ib2, O(1), O(2), O(3), dE
        if (i1.ne.ib1 .or. i2.ne.ib2) then
           print *, 'ERROR', i1, ib1, i2, ib2
        endif
        !if (i1.GE.nemin .AND. i1.LT.nemin+nbands .AND. i2.GE.nemin .AND. i2.LT.nemin+nbands) then
        if (i1.LE.NE .and. i2.LE.NE) then
           DO ixyz=1,3
              Vel(i1,i2,ixyz) = O(ixyz)
              if (i2.NE.i1) then
                 Vel(i2,i1,ixyz) = conjg(O(ixyz))
              endif
           ENDDO
        endif
     enddo
  enddo
  if (nb_max.GT.NE) nb_max = NE
!9010 FORMAT(/,2X,' KP:',I6,' NEMIN NEMAX : ',2I5, ' dE:',2f5.2,' K:',a10 /)
9010 FORMAT(/,2X,4X,I6,15X,2I5,4X,2f5.2,3X,a10 /)
9040 FORMAT(3X,2I4,6E13.6,F13.8)
END SUBROUTINE ReadVelocityK

SUBROUTINE ReadEnergiesK(Ek, wgh, NE, ikp, nkpt, nume)
  IMPLICIT NONE
!!!! Reading band energies from energy file
  REAL*8, intent(out) :: wgh(nkpt), Ek(nume)
  INTEGER, intent(out):: NE
  INTEGER, intent(in) :: ikp, nkpt, nume
  ! locals
  INTEGER :: is, ios, itape, ii, NUM
  REAL*8  :: S, T, Z, WG, E1
  CHARACTER*10 :: KNAME
  INTEGER :: N
  itape=59
  READ(itape,'(3e19.12,a10,2i6,e19.12)',IOSTAT=ios) S,T,Z,KNAME,N,NE,WG
  wgh(ikp) = WG
  DO ii=1,NE
     READ(itape,*) NUM, E1
     Ek(NUM)=E1
  ENDDO
END SUBROUTINE ReadEnergiesK

SUBROUTINE ReadDMFT_TransK_Outside(nsymop, nbands, nemin, numk, wkii, wkis, fUdmft, fh_p, pform, ikp)
!!! Start reading DMFT transformation
  IMPLICIT NONE
  INTEGER, intent(out)   :: nsymop, nbands, nemin
  INTEGER, intent(inout) :: numk, wkii, wkis
  CHARACTER*100, intent(in) :: fUdmft
  INTEGER, intent(in)    :: fh_p, ikp
  LOGICAL, intent(in)    :: pform
  ! locals
  CHARACTER*10 :: skii
  INTEGER :: tnorbitals, iikp, tmaxdim2
  if (wkii==numk) then ! There are many udmfile's because of parallel execution! Need to read one after the other.
     wkii=0
     wkis = wkis + 1
     WRITE(skii,fmt='(I2)') wkis
     close(fh_p)
     if (pform) then
        open(fh_p, file=TRIM(fUdmft)//ADJUSTL(TRIM(skii)), status='old', form='formatted')
        READ(fh_p,*) numk, nsymop, tnorbitals
     else
        open(fh_p, file=TRIM(fUdmft)//ADJUSTL(TRIM(skii)), status='old', form='unformatted')
        READ(fh_p) numk, nsymop, tnorbitals
     endif
  endif

  if (pform) then
     READ(fh_p,*) iikp, nbands, tmaxdim2, tnorbitals, nemin
  else
     READ(fh_p) iikp, nbands, tmaxdim2, tnorbitals, nemin
  endif

  if (iikp.NE.ikp) print *, 'ERROR: ikp and iikp=', ikp, iikp

  wkii = wkii + 1
END SUBROUTINE ReadDMFT_TransK_Outside

SUBROUTINE ReadDMFT_TransK_Inside(DMFTU, fh_p, pform, isym, nindo, nbands, norbitals, maxdim2)
  ! Continuing reading DMFT transformation
  IMPLICIT NONE
  COMPLEX*16, intent(out) :: DMFTU(nbands,maxdim2,norbitals)
  INTEGER, intent(in)     :: fh_p, isym, nindo(norbitals), nbands, norbitals, maxdim2
  LOGICAL, intent(in)     :: pform
  ! local
  INTEGER :: iisym, iorb, ind, i
  if (pform) then
     READ(fh_p,*) iisym
  else
     READ(fh_p) iisym
  endif
  if (iisym.NE.isym) print *, 'ERROR: isym and iisym=', isym, iisym
  DO iorb=1,norbitals
     DO ind=1,nindo(iorb)
        if (pform) then
           !READ(fh_p,*) (DMFTU(i,ind,iorb),i=1,nbands)
           DO i=1,nbands
              READ(fh_p,'(f16.9,1x,f16.9,3x)',advance='no') DMFTU(i,ind,iorb)
           ENDDO
           READ(fh_p,*)
        else
           READ(fh_p) (DMFTU(i,ind,iorb),i=1,nbands)
        endif
     ENDDO
  ENDDO
  !DO iorb=1,norbitals
  !   DO ind=1,nindo(iorb)
  !      DO i=1,nbands
  !         WRITE(6,'(I2,1x,I2,1x,I2,1x,I3,1x,2f14.9)') iisym,iorb,ind,i,DMFTU(i,ind,iorb)
  !      ENDDO
  !   ENDDO
  !ENDDO
END SUBROUTINE ReadDMFT_TransK_Inside

SUBROUTINE Read_The_Rest_Basic_Arrays(fhb, nindo, cixdim, nl, ll, iorbital, cix, nind, csize, iSx, Sigind, EF, VOL, norbitals, ncix, natom, maxdim, maxdim2)
  IMPLICIT NONE
  INTEGER, intent(in)  :: fhb
  INTEGER, intent(out) :: nindo(norbitals), cixdim(ncix), nl(natom)
  INTEGER, intent(out) :: ll(natom,4), iorbital(natom,4), cix(natom,4), nind(natom,4)
  INTEGER, intent(out) :: csize(ncix), iSx(maxdim2,norbitals), Sigind(maxdim,maxdim,ncix)
  REAL*8, intent(out)  :: EF, VOL
  INTEGER, intent(in)  :: norbitals, ncix, natom, maxdim, maxdim2
  ! locals
  character*100 :: STR
  INTEGER :: iorb, icix, icase, lcase, ip, iq, ip1
  !!! ReadBasicArrays
  READ(fhb, *) STR ! 'nindo'
  READ(fhb, *) (nindo(iorb),iorb=1,norbitals)
  READ(fhb, *) STR ! 'cidim'
  READ(fhb, *) (cixdim(icix),icix=1,ncix)
  READ(fhb, *) STR ! 'nl'
  READ(fhb, *) (nl(icase), icase=1,natom)
  READ(fhb, *) STR ! 'll, iorbital, cix, nind'
  do icase=1,natom
     do lcase=1,nl(icase)
        READ(fhb, *) ll(icase,lcase)
        READ(fhb, *) iorbital(icase,lcase)
        READ(fhb, *) cix(icase,lcase)
        READ(fhb, *) nind(icase,lcase)          ! nind
     enddo
  enddo

  READ(fhb, *) STR ! csize
  READ(fhb, *) (csize(icix), icix=1,ncix)
  READ(fhb, *) STR ! 'iSx'
  do iorb=1,norbitals
     do ip1=1,nindo(iorb)
        READ(fhb, *) iSx(ip1,iorb)
     enddo
  enddo
  READ(fhb, *) STR ! 'Sigind'
  DO icix=1,ncix
     DO ip=1,cixdim(icix)
        DO iq=1,cixdim(icix)
           READ(fhb, *) Sigind(ip,iq,icix)
        ENDDO
     ENDDO
  ENDDO
  READ(fhb, *) STR ! EF
  READ(fhb, *) EF, VOL
END SUBROUTINE Read_The_Rest_Basic_Arrays

SUBROUTINE Transform_To_Reducible(Ve, isym, Vel, opimat, NE, nord)
  IMPLICIT NONE
  COMPLEX*16, intent(out) :: Ve(NE,NE,3)
  INTEGER, intent(in)     :: isym
  COMPLEX*16, intent(in)  :: Vel(NE,NE,3)
  REAL*8, intent(in)      :: opimat(3,3,nord)
  INTEGER, intent(in)     :: NE, nord
  ! locals
  INTEGER    :: ip, iq, i, ii
  COMPLEX*16 :: osm
  ! Changing matrix elements to general k-point
  DO ip=1,NE
     DO iq=1,NE
        do i=1,3
           osm=0
           do ii=1,3
              osm = osm + Vel(ip,iq,ii)*opimat(ii,i,isym)
           end do
           Ve(ip,iq,i) = osm
        end do
     ENDDO
  ENDDO
END SUBROUTINE Transform_To_Reducible

SUBROUTINE Opt_Inside_Simple(copt, Cv, Dv, zekw, k, l, amesh, iw_p, iw_m, wdelta, wgamma, EF, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, nb_min, nb_max, InterbandOnly)
  IMPLICIT NONE
  COMPLEX*16, intent(out):: copt(Ndirection)
  COMPLEX*16, intent(in) :: Cv(NE,NE,Ndirection), Dv(NE,NE,Ndirection), zekw(NE,2*Nd+1)
  INTEGER, intent(in) :: k, l, amesh(2*Nd+1), iw_p, iw_m, Ndirection
  REAL*8, intent(in)  :: wdelta, wgamma, EF
  INTEGER, intent(in) :: nbands, Nd, bandind(NE), nume, NE, nemin, nb_min, nb_max
  LOGICAL, intent(in) :: InterbandOnly, Qinside(NE)
  ! locals
  REAL*8  :: dx, wdx
  INTEGER :: p, q, dir
  COMPLEX*16  :: eps_p, eps_m, womega, IMAG
  !------------------------------------------------------------------
  DATA IMAG/(0.0D0,1.0D0)/

  eps_p = amesh(iw_p)*wdelta    ! eps
  eps_m = amesh(iw_m)*wdelta    ! eps-omega
  womega = amesh(k+Nd+1)*wdelta ! omega

  dx = womega/k
  !dx=(amesh(Nd+2)-amesh(Nd+1))*wdelta
  wdx=dx
  if (l.EQ.1 .or. l.EQ.k+1) wdx=0.5*dx ! Trapezoid rule integration: the first and last interval has 1/2 of the integration weight
  copt = 0
  DO dir=1,Ndirection
     DO p=nb_min,nb_max
        DO q=nb_min,nb_max
           if (InterbandOnly .and. bandind(p).eq.bandind(q)) CYCLE
           if (.NOT.Qinside(p) .AND. .NOT.Qinside(q) ) CYCLE  !! both p and q are outside the window -> non-interacting problem computed outside
           copt(dir) = copt(dir) + wdx*Cv(p,q,dir)/(     (eps_p+EF-zekw(p,iw_p)+wgamma*IMAG)*(eps_m+EF-zekw(q,iw_m)+wgamma*IMAG))
           copt(dir) = copt(dir) - wdx*Dv(p,q,dir)/(conjg(eps_p+EF-zekw(p,iw_p)+wgamma*IMAG)*(eps_m+EF-zekw(q,iw_m)+wgamma*IMAG))
        ENDDO
     ENDDO
  ENDDO
END SUBROUTINE Opt_Inside_Simple

SUBROUTINE Opt_Inside(copt, Cv, Dv, zekw, k, l, amesh, iw_p, iw_m, wdelta, wgamma, EF, lg_q, lg_p, nbands, Nd, bandind, Qinside, nume, Ndirection, NE, nemin, nb_min, nb_max, InterbandOnly)
  ! Function computes:
  !          Cv(ip,iq)*wes - Dv(ip,iq)*wos
  !          where wes and wos are appropraite weights obtained by integrating over frequency
  !          in a small interval of frequency. Frequency points are given by amesh(:)*wdelta
  IMPLICIT NONE
  COMPLEX*16, intent(out)   :: copt(Ndirection)
  COMPLEX*16, intent(in)    :: Cv(NE,NE,Ndirection), Dv(NE,NE,Ndirection), zekw(NE,2*Nd+1)
  INTEGER, intent(in)       :: k, l, amesh(2*Nd+1), iw_p, iw_m
  REAL*8, intent(in)        :: wdelta, wgamma, EF
  COMPLEX*16, intent(inout) :: lg_q(NE), lg_p(NE)
  INTEGER, intent(in)       :: nbands, Nd, bandind(NE), NE, nume, Ndirection, nb_min, nb_max, nemin
  LOGICAL, intent(in)       :: InterbandOnly, Qinside(NE)
  ! locals
  REAL*8  :: dx, wdx, a, b, womega, eps_p, eps_m, big
  INTEGER :: p, q, dir
  COMPLEX*16 :: denom_es, denom_os, wes, wos
  COMPLEX*16 :: IMAG
  !------------------------------------------------------------------
  DATA IMAG/(0.0D0,1.0D0)/

  big = 1e6

  eps_p = amesh(iw_p)*wdelta    ! eps
  eps_m = amesh(iw_m)*wdelta    ! eps-omega
  womega = amesh(k+Nd+1)*wdelta ! omega

  if (l.EQ.1 .or. iw_p.EQ.1) then        ! The first interval to integrate is smaller
     a=amesh(iw_p)*wdelta
  else
     a=0.5*(amesh(iw_p-1)+amesh(iw_p))*wdelta
  endif
  if (l.eq.k+1 .or. iw_p.EQ.2*Nd+1) then ! The last interval to integrate is smaller
     b=amesh(iw_p)*wdelta
  else                    ! The rest of the intervals are the same
     b=0.5*(amesh(iw_p)+amesh(iw_p+1))*wdelta
  endif

  DO p=nb_min,nb_max
     lg_q(p) = log(b+EF-womega-zekw(p,iw_m)+wgamma*IMAG) - log(a+EF-womega-zekw(p,iw_m)+wgamma*IMAG)
     lg_p(p) = log(b+EF-       zekw(p,iw_p)+wgamma*IMAG) - log(a+EF-       zekw(p,iw_p)+wgamma*IMAG)
  ENDDO

  copt=0
  DO dir=1,Ndirection
     DO p=nb_min,nb_max
        DO q=nb_min,nb_max
           if (InterbandOnly .and. bandind(p).eq.bandind(q)) CYCLE
           if (.NOT.Qinside(p) .AND. .NOT.Qinside(q) ) CYCLE  !! both p and q are outside the window -> non-interacting problem computed outside
           denom_es = womega -       zekw(p,iw_p)  + zekw(q,iw_m)
           denom_os = womega - conjg(zekw(p,iw_p)) + zekw(q,iw_m) -2*wgamma*IMAG
           wes = (lg_q(q)-      lg_p(p)) /denom_es
           wos = (lg_q(q)-conjg(lg_p(p)))/denom_os
           if (abs(wes).gt.big) then
              print *, 'Warning: the weightes is very big at ', womega, eps_p, p, q, wes
              wes=0
           endif
           if (abs(wos).gt.big) then
              print *, 'Warning: the weightos is very big at ', womega, eps_p, p, q, wos
              wos=0
           endif
           copt(dir) = copt(dir) + Cv(p,q,dir)*wes - Dv(p,q,dir)*wos
        ENDDO
     ENDDO
  ENDDO
  RETURN
END SUBROUTINE Opt_Inside

SUBROUTINE FindEigsys(zek, Al, Ar, ix, wsigma, gij, womega, omega, STrans, Ek, nom, ncix, nbands, nemin, maxsize, nume, csize, sigma)
  IMPLICIT NONE
  COMPLEX*16, intent(out)   :: zek(nbands), Al(nbands,nbands), Ar(nbands, nbands)
  INTEGER, intent(inout)    :: ix
  COMPLEX*16, intent(inout) :: wsigma(maxsize, ncix)
  COMPLEX*16, intent(inout) :: gij(nbands,nbands)
  REAL*8, intent(in)     :: womega
  REAL*8, intent(in)     :: omega(nom)
  COMPLEX*16, intent(in) :: STrans(maxsize,ncix,nbands,nbands)
  REAL*8, intent(in)     :: Ek(nume)
  INTEGER, intent(in)    :: nom, ncix, nbands, nemin, maxsize, nume, csize(ncix)
  COMPLEX*16, intent(in) :: sigma(nom, maxsize, ncix)
  ! locals
  INTEGER :: icix, ind, i

  gij=0
  DO i=1,nbands
     gij(i,i) = Ek(i+nemin-1)  !-----  g^-1 of the LDA part in band representation ----!
  ENDDO
  CALL findNext(womega, ix, omega, nom)
  do icix=1,ncix
     do ind=1,csize(icix)
        CALL interp(wsigma(ind,icix), sigma(:,ind,icix), womega, ix, omega, nom)
     enddo
  enddo
  CALL AddSigma_optimized2(gij, wsigma, STrans, csize, 1, nbands, ncix, maxsize)
  CALL eigsys(gij, zek, Al, Ar, nbands)
END SUBROUTINE FindEigsys

SUBROUTINE cmp_C_and_D(Cv, Dv, Al_p, Ar_p, Al_m, Ar_m, Ve, alphaV, Ndirection, L1, L2, L3, L4, Ltmp, Qinside, ndim, NE, nemin, nemax, nb_min, nb_max)
  !--------------------------------------
  IMPLICIT NONE
  !---------- Passed variables ----------
  COMPLEX*16, intent(out)   :: Cv(NE,NE,Ndirection), Dv(NE,NE,Ndirection)
  COMPLEX*16, intent(in)    :: Al_p(ndim,ndim), Ar_p(ndim,ndim), Al_m(ndim,ndim), Ar_m(ndim,ndim)
  COMPLEX*16, intent(in)    :: Ve(NE,NE,3)
  REAL*8, intent(in)        :: alphaV(3,3,Ndirection)
  COMPLEX*16, intent(inout) :: L1(NE,NE,3), L2(NE,NE,3), L3(NE,NE,3), L4(NE,NE,3), Ltmp(ndim,ndim)
  INTEGER, intent(in)       :: ndim, Ndirection, NE, nb_min, nb_max, nemin, nemax
  LOGICAL, intent(in)       :: Qinside(NE)
  !---------- Local variables -----------
  COMPLEX*16 :: Al_p_c(ndim,ndim), Ar_p_c(ndim,ndim)
  COMPLEX*16 :: c1tmp, c2tmp
  INTEGER    :: p,q, i,j,l, dir

  ! Calculation of matrices C and D that are multiplied to the double denominator
  ! C is for equal sign and D for opposite sign
  Al_p_c = conjg(transpose(Al_p))
  Ar_p_c = conjg(transpose(Ar_p))

  DO i=1,3
     ! inside the small window
     Ltmp = matmul(Al_m, Ve(nemin:nemax,nemin:nemax,i))
     L1(nemin:nemax,   nemin:nemax,   i) = matmul(Ltmp, Ar_p)
     L3(nemin:nemax,   nemin:nemax,   i) = matmul(Ltmp, Al_p_c)

     ! cross terms, where one index is inside window, and one is outside
     if (nb_min.LT.nemin) then
        L1(nb_min:nemin-1,nemin:nemax,   i) = matmul(Ve(nb_min:nemin-1,nemin:nemax,i),Ar_p)
        L1(nemin:nemax,   nb_min:nemin-1,i) = matmul(Al_m, Ve(nemin:nemax,nb_min:nemin-1,i))
        L3(nb_min:nemin-1,nemin:nemax,   i) = matmul(Ve(nb_min:nemin-1,nemin:nemax,i), Al_p_c)
        L3(nemin:nemax,   nb_min:nemin-1,i) = L1(nemin:nemax,nb_min:nemin-1,i)
     endif
     if (nb_max.GT.nemax) then
        L1(nemax+1:nb_max,nemin:nemax,   i) = matmul(Ve(nemax+1:nb_max,nemin:nemax,i),Ar_p)
        L1(nemin:nemax,   nemax+1:nb_max,i) = matmul(Al_m, Ve(nemin:nemax,nemax+1:nb_max,i))
        L3(nemax+1:nb_max,nemin:nemax,   i) = matmul(Ve(nemax+1:nb_max,nemin:nemax,i), Al_p_c)
        L3(nemin:nemax,nemax+1:nb_max,   i) = L1(nemin:nemax,nemax+1:nb_max,i)
     endif

     ! inside the small window
     Ltmp = matmul(Ve(nemin:nemax,nemin:nemax,i), Ar_m)
     L2(nemin:nemax,   nemin:nemax,   i) = matmul(Al_p, Ltmp)
     L4(nemin:nemax,nemin:nemax,i) = matmul(Ar_p_c, Ltmp)

     ! cross terms, where one index is inside window, and one is outside
     if (nb_min.LT.nemin) then
        L2(nb_min:nemin-1,nemin:nemax,   i) = matmul(Ve(nb_min:nemin-1,nemin:nemax,i),Ar_m)
        L2(nemin:nemax,   nb_min:nemin-1,i) = matmul(Al_p, Ve(nemin:nemax,nb_min:nemin-1,i))
        L4(nb_min:nemin-1,nemin:nemax,   i) = L2(nb_min:nemin-1,nemin:nemax,   i)
        L4(nemin:nemax,   nb_min:nemin-1,i) = matmul(Al_p_c, Ve(nemin:nemax,nb_min:nemin-1,i))
     endif
     if (nb_max.GT.nemax) then
        L2(nemax+1:nb_max,nemin:nemax,   i) = matmul(Ve(nemax+1:nb_max,nemin:nemax,i),Ar_m)
        L2(nemin:nemax,   nemax+1:nb_max,i) = matmul(Al_p, Ve(nemin:nemax,nemax+1:nb_max,i))
        L4(nemax+1:nb_max,nemin:nemax,   i) = L2(nemax+1:nb_max,nemin:nemax,   i)
        L4(nemin:nemax,   nemax+1:nb_max,i) = matmul(Al_p_c, Ve(nemin:nemax,nemax+1:nb_max,i))
     endif
  ENDDO

  DO p=nb_min,nb_max
     DO q=nb_min,nb_max
        ! When both p and q are outside, we do not need to update Cv and Dv, because they
        ! are equal to v^2.
        if (.NOT.Qinside(p) .AND. .NOT.Qinside(q) ) CYCLE
        DO dir=1,Ndirection
           c1tmp=0
           c2tmp=0
           DO j=1,3
              DO l=1,3
                 c1tmp = c1tmp + L2(p,q,j)*L1(q,p,l)*alphaV(l,j,dir)
                 c2tmp = c2tmp + L4(p,q,j)*L3(q,p,l)*alphaV(l,j,dir)
              ENDDO
           ENDDO
           Cv(p,q,dir) = c1tmp
           Dv(p,q,dir) = c2tmp
        ENDDO
     ENDDO
  ENDDO
END SUBROUTINE cmp_C_and_D

SUBROUTINE cmp_C_and_D_old(Cv, Dv, Al_p, Ar_p, Al_m, Ar_m, Ve, alphaV, Ndirection, L1, L2, L3, L4, Ltmp, ndim)
  !--------------------------------------
  IMPLICIT NONE
  !---------- Passed variables ----------
  COMPLEX*16, intent(out)   :: Cv(ndim,ndim,Ndirection), Dv(ndim,ndim,Ndirection)
  COMPLEX*16, intent(in)    :: Al_p(ndim,ndim), Ar_p(ndim,ndim), Al_m(ndim,ndim), Ar_m(ndim,ndim)
  COMPLEX*16, intent(in)    :: Ve(ndim,ndim,3)
  REAL*8, intent(in)        :: alphaV(3,3,Ndirection)
  COMPLEX*16, intent(inout) :: L1(ndim,ndim,3), L2(ndim,ndim,3), L3(ndim,ndim,3), L4(ndim,ndim,3), Ltmp(ndim,ndim)
  INTEGER, intent(in)       :: ndim, Ndirection
  !---------- Local variables -----------
  COMPLEX*16 :: c1tmp, c2tmp
  INTEGER    :: p,q, i,j,l, dir

  ! Calculation of matrices C and D that are multiplied to the double denominator
  ! C is for equal sign and D for opposite sign
  DO i=1,3
     Ltmp = matmul(Al_m, Ve(:,:,i))
     L1(:,:,i) = matmul(Ltmp, Ar_p)
     L3(:,:,i) = matmul(Ltmp, conjg(transpose(Al_p)))
     Ltmp = matmul(Ve(:,:,i), Ar_m)
     L2(:,:,i) = matmul(Al_p, Ltmp)
     L4(:,:,i) = matmul(conjg(transpose(Ar_p)), Ltmp)
  ENDDO
  Cv = 0
  Dv = 0
  DO p=1,ndim
     DO q=1,ndim
        DO dir=1,Ndirection
           c1tmp=0
           c2tmp=0
           DO j=1,3
              DO l=1,3
                 c1tmp = c1tmp + L2(p,q,j)*L1(q,p,l)*alphaV(l,j,dir)
                 c2tmp = c2tmp + L4(p,q,j)*L3(q,p,l)*alphaV(l,j,dir)
              ENDDO
           ENDDO
           Cv(p,q,dir) = Cv(p,q,dir) + c1tmp
           Dv(p,q,dir) = Dv(p,q,dir) + c2tmp
        ENDDO
     ENDDO
  ENDDO

END SUBROUTINE cmp_C_and_D_old

SUBROUTINE Opt_Noncorrelated(copt, Cv0, Ek, womega, EF, wgamma,nb_min, nb_max, nemin, nemax, nume, Qinside, NE, Ndirection, InterbandOnly)
!!! Add non-correlated bands !!!
  IMPLICIT NONE
!!! For non-correlated
  COMPLEX*16, intent(out) :: copt(Ndirection)
  COMPLEX*16, intent(in) :: Cv0(NE,NE,Ndirection)
  REAL*8, intent(in)     :: Ek(nume)
  REAL*8, intent(in)  :: womega, EF
  REAL*8, intent(in)  :: wgamma
  INTEGER, intent(in) :: nb_min, nb_max, nemin, nemax, nume, NE, Ndirection
  LOGICAL, intent(in) :: InterbandOnly, Qinside(NE)
  ! local variables
  REAL*8     :: aw, bw
  INTEGER    :: ip, iq, dir
  COMPLEX*16 :: denom_es, denom_os, wes, wos
  COMPLEX*16, ALLOCATABLE :: lg_p(:), lg_q(:)
  COMPLEX*16, PARAMETER  :: IMAG = (0.0D0,1.0D0)

  ALLOCATE( lg_p(nb_max), lg_q(nb_max) )

  aw=0.0
  bw=womega
  DO ip=nb_min,nb_max
     if (Qinside(ip)) CYCLE ! These bands were taken into account in DMFT part
     lg_q(ip) = log(bw+EF-womega-Ek(ip)+wgamma*IMAG) - log(aw+EF-womega-Ek(ip)+wgamma*IMAG)
     lg_p(ip) = log(bw+EF-       Ek(ip)+wgamma*IMAG) - log(aw+EF-       Ek(ip)+wgamma*IMAG)
  ENDDO
  copt=0
  DO ip=nb_min,nb_max

     !if (InterbandOnly .OR. (ip.GE.nemin .and.ip.LE.nemax)) CYCLE ! These bands were taken into account in DMFT part
     DO iq=nb_min,nb_max

        if ( Qinside(ip) .or. Qinside(iq) ) CYCLE ! These bands were taken into account in DMFT part

        denom_es = womega - Ek(ip) + Ek(iq)
        denom_os = womega - Ek(ip) + Ek(iq) -2*wgamma*IMAG
        wes = (lg_q(iq)-      lg_p(ip)) /denom_es
        wos = (lg_q(iq)-conjg(lg_p(ip)))/denom_os
        DO dir=1,Ndirection
           copt(dir) = copt(dir) + Cv0(ip,iq,dir)*(wes - wos)
        ENDDO
     ENDDO
  ENDDO
  DEALLOCATE( lg_p, lg_q )
END SUBROUTINE Opt_Noncorrelated

program read_optm
  USE com_mpi, ONLY: start_MPI, stop_MPI
  IMPLICIT NONE
  ! functions
  INTEGER :: CountSelfenergy
  ! variables
  CHARACTER*100 :: case, fmommat, fsymop, fUdmft, fenergy, fbasicArrays, fsigname
  CHARACTER*2   :: updn, dnup, so
  INTEGER :: fhb, fhi, fh_m
  CHARACTER*100 :: STR
  INTEGER :: nom, nk, iso, nat, norbitals, ncix, natom, nkpt, nmat, nume, lmax2, maxdim2, maxdim, maxsize
  LOGICAL :: Qsym, Qsimple, Qcomplex, InterbandOnly, ProjectToCorrelated
  REAL*8  :: gamma, gammac, ommax, delta, dwindow, Temperature
  INTEGER :: Nd, Nd0, i, j, l
  CHARACTER*5 :: adum
  INTEGER :: Ndirection
  LOGICAL :: file_exists
  REAL*8, ALLOCATABLE :: alphaV(:,:,:)

  CALL start_MPI()

  Qsimple = .False.
  updn=''
  fhi = 995
  open(fhi, file='dmftopt.in', status='old', form='formatted')
  READ(fhi, *) Temperature
  READ(fhi, *) case
  READ(fhi, *) gamma
  READ(fhi, *) gammac
  READ(fhi, *) ommax
  READ(fhi, *) delta
  READ(fhi, *) Nd0
  READ(fhi, *) Qsym
  READ(fhi, *) InterbandOnly
  READ(fhi, *) dwindow
  READ(fhi, *) Ndirection
  print *, 'Ndirection=', Ndirection
  ALLOCATE( alphaV(3,3,Ndirection) )
  do l=1,Ndirection
     READ(fhi, *) alphaV(1,1,l), alphaV(1,2,l), alphaV(1,3,l)
     READ(fhi, *) alphaV(2,1,l), alphaV(2,2,l), alphaV(2,3,l)
     READ(fhi, *) alphaV(3,1,l), alphaV(3,2,l), alphaV(3,3,l)
  enddo
  close(fhi)

  !if (.NOT.(updn.EQ.'up' .OR.updn.EQ.'dn')) updn=''
  !!
  !dnup='dn'
  !IF (updn.EQ.'dn') THEN
  !   dnup = 'up'
  !ELSEIF (updn.EQ.'up') THEN
  !   dnup = 'dn'
  !ELSE
  !   updn = ''
  !ENDIF
  !!

  !!! Counting the number of frequency points in self-energy
  fsigname = 'sig.inp'//TRIM(updn)//'1'
  print *, 'fsigname=', fsigname
  open(81, file=fsigname, status='old', form='formatted')
  nom = CountSelfenergy(81, 1)-1 !--- How many frequency point exists in the input file? ---!
  close(81)

  fbasicArrays = 'BasicArrays.dat'
  !!! Reading the first part of the header file to get dimensions
  fhb=996
  open(fhb,file=fbasicArrays,status='old', form='formatted')
  READ(fhb, *) STR ! 'nat, iso, norbitals, ncix, natom'
  READ(fhb, *) nat, iso, norbitals, ncix, natom
  READ(fhb, *) STR ! nkpt, nmat, nume
  READ(fhb, *) nkpt, nmat, nume
  READ(fhb, *) STR ! Qcomplex
  READ(fhb, *) Qcomplex
  READ(fhb, *) STR ! 'lmax2, maxdim2, maxdim, maxsize'
  READ(fhb, *) lmax2, maxdim2, maxdim, maxsize

  so = ''
  IF (iso.EQ.2) so='so'
  fenergy = TRIM(case)//'.energy'//TRIM(so)//TRIM(updn)
  fmommat = TRIM(case)//'.mommat'//TRIM(updn)
  fsymop = TRIM(case)//'.symop'
  fUdmft = 'Udmft'//TRIM(updn)//'.'

  fh_m = 99  ! file for velocity matrix elements

  INQUIRE(file=fmommat, exist=file_exists)
  if (.not. file_exists) then
     fmommat = TRIM(case)//'.mommat2'//TRIM(updn)
     INQUIRE(file=fmommat, exist=file_exists)
     if (.not. file_exists) WRITE(*,*) 'ERROR File not found: Expecting that ', fmommat, 'exists.'
  endif
  open(fh_m, file=fmommat, status='old')

  if (Temperature.NE.0.0) then
     ! We will use linear mesh only. We need to set very large Nd, so that
     ! cutoff >=4*T, and ommax/(2*delta)=Nd
     if (Nd0 .LT. int(2*Temperature/Delta+1.)) Nd0=int(2*Temperature/Delta+1.)
     if (ommax.LT.2*Nd0*delta) then
        ommax=2*Nd0*delta
     else
        Nd0 = ommax/(2*delta)
     endif
  endif

  WRITE(*,'(A,A)') 'case=', TRIM(case)
  WRITE(*,'(A,f10.7,1x,A,f10.7)') 'gamma=', gamma, 'gammac=', gammac
  WRITE(*,'(A,f10.3,1x,A,f14.10,1x,A,I4)') 'ommax=', ommax, 'delta=', delta, 'Nd=', Nd0
  WRITE(*,'(A,L)') 'Qsym=', Qsym
  WRITE(*, '(A)') 'alphaV='
  WRITE(*,*) (((alphaV(j,i,l), i=1,3), j=1,3),l=1,Ndirection)
  WRITE(*,'(A,A)') 'fmommat=', fmommat
  WRITE(*,'(A,A)') 'fenergy=', fenergy
  WRITE(*,'(A,A)') 'fUdmft=', fUdmft
  WRITE(*,'(A,A)') 'fsymop=', fsymop
  WRITE(*,'(A,L)') 'InterbandOnly=', InterbandOnly
  WRITE(*,'(A,f10.7)') 'dwindow=', dwindow
  WRITE(*,'(A,f10.7)') 'Temperature=', Temperature

  Nd = 2*Nd0
  CALL Cmp_Optics(fenergy, fUdmft, fsymop, fh_m, fhb, ommax, Temperature, delta, Nd, gammac, gamma, alphaV, Ndirection, Qsym, Qsimple, nom, nat, iso, norbitals, ncix, natom, nkpt, nmat, nume, Qcomplex, lmax2, maxdim2, maxdim, maxsize, dwindow, InterbandOnly)

  close(fh_m)
  DEALLOCATE( alphaV )
  CALL stop_MPI()

end program read_optm

REAL*8 Function ferm(x)
  IMPLICIT NONE
  REAL*8, intent(in):: x
  if (x.gt.100) then
     ferm=0.0
     return
  endif
  if (x.lt.-100) then
     ferm=1.0
     return
  endif
  ferm = 1./(1.+exp(x))
  return
end function ferm


text:
! cable_air.f90
!
! Source file containing air routines for CABLE
!
! Development by Ying-Ping Wang, Eva Kowalczyk, Ray Leuning
! Gab Abramowitz, Martin Dix, Harvey Davies, Mike Raupach,
!
! bugs to bernard.pak@csiro.au
!
! This file contains modules:
!   air_module,
! The subroutines included are:
!   define_air,
!
! Most user-defined types (e.g. met%tk) are defined in cable_types module
! in cable_variables.f90

MODULE cable_air
  USE cable_physical_constants
  USE cable_dimensions, ONLY: mp_patch, r_1
  USE cable_types
  IMPLICIT NONE
  PRIVATE
  PUBLIC define_air
CONTAINS
  !-----------------------------------------------------------------------
  SUBROUTINE define_air(met,air)
    TYPE (air_type), INTENT(INOUT) :: air ! air_type variables
    TYPE (met_type), INTENT(IN)  :: met ! meteorological variables
    REAL(r_1), DIMENSION(mp_patch) :: es ! sat vapour pressure (mb)
    ! Calculate saturation vapour pressure
    es = tetena * EXP(tetenb * (met%tvair-tfrz)/(tetenc + (met%tvair-tfrz)))
    ! Calculate conversion factor from from m/s to mol/m2/s
    air%cmolar = met%pmb * 100.0 / (rgas * (met%tvair))
    ! Calculate dry air density:
    air%rho = MIN(1.3,rmair * air%cmolar)
    ! molar volume (m^3/mol)
    air%volm = rgas * (met%tvair) / (100.0 * met%pmb)
    ! latent heat for water (j/kg)
    ! air%rlam = (2501.0 - 2.38 * (met%tvair- tfrz)) * 1000.0
	air%rlam  = 1.91846e6*((met%tvair)/((met%tvair)-33.91))**2  ! Henderson-Sellers, QJRMS, 1984
    !air%rlam=2.5104e6
    ! saturation specific humidity
    air%qsat = (rmh2o / rmair) * es / met%pmb
    ! d(qsat)/dT ((kg/kg)/K)
    air%epsi = (air%rlam / capp) * (rmh2o / rmair) * es * tetenb * tetenc / &
         &  (tetenc + (met%tvair - tfrz)) ** 2 / met%pmb
    ! air kinematic viscosity (m^2/s)
    air%visc = 1e-5 * MAX(1.0, 1.35 + 0.0092 * (met%tvair - tfrz))
    ! psychrometric constant
    air%psyc = met%pmb * 100.0 * capp * rmair / air%rlam / rmh2o
    ! d(es)/dT (mb/K)
    air%dsatdk = 100.0*(tetena*tetenb*tetenc)/((met%tvair-tfrz)+tetenc)**2 &
         * EXP(tetenb*(met%tvair-tfrz)/((met%tvair-tfrz) + tetenc))
    !      air%dsatdk = (610.078 * 17.27 * 237.3) / ((met%tvair-tfrz)+237.2)** 2 * &
    !                 & EXP(17.27 * (met%tvair-tfrz) / ((met%tvair-tfrz) + 237.3))
  END SUBROUTINE define_air
END MODULE cable_air

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


module hamiltonian

  use precision
  use mpi
  use globals
  use fft
  use misc
  use potential

  implicit none

contains

  function f_hamiltonian(wavefn, nl, nu) result(hamiltonian)

    !
    ! This function computes the hamiltonian of a wave function
    !

    integer(idp), dimension(2), intent(in) :: nl, nu
    complex(cdp), dimension(nl(1):nu(1), nl(2):nu(2)), intent(in) ::&
         & wavefn
    complex(cdp), dimension(nl(1):nu(1), nl(2):nu(2)) :: hamiltonian

    complex(cdp), dimension(:, :), allocatable, save :: del2, pot,&
         & delwavefnx, delwavefny, delwavefn, delx, dely, delxy,&
         & vecpotfn2
    real(fdp), dimension(:, :), allocatable, save :: vecpotfnx,&
         & vecpotfny, vecpotfn
    integer(idp), dimension(dim), save :: nlold = 0, nuold = 0
    logical, save :: defvecpot = .false., defpot = .false., defdel =&
         & .false.

    !
    ! If the space size changes, deallocate all the stored function
    !

    if (nlold(1) /= nl(1) .or. nlold(2) /= nl(2) .or. nuold(1) /=&
         & nu(1) .or. nuold(2) /= nu(2)) then

       if (defdel) then

          deallocate(del2)

          defdel = .false.

       end if

       if (defvecpot) then

          select case (parameters%gauge)

          case('(-y,0,0)')
             deallocate(delx, vecpotfny, vecpotfn2, delwavefnx)

          case('(-y,x,0)')

             deallocate(delx, dely, vecpotfnx,  vecpotfny, delwavefnx&
                  &, delwavefny, vecpotfn2)

          case('(x-y,x-y,0)', 'fft')
             deallocate(delxy, vecpotfn,  delwavefn, vecpotfn2)

          end select

          defvecpot = .false.

       end if

       if (defpot) then

          deallocate(pot)

          defpot = .false.

       end if

       nlold = nl
       nuold = nu

    end if

    !
    ! If the required function are not allocated and computed,
    ! allocate them and compute them
    !

    if (.not. defdel) then

       allocate(del2(nl(1):nu(1), nl(2):nu(2)))

       call s_momentumfac(del2, nl, nu, (/ 1.0d0, 1.0d0 /)&
            &/dble(product(parameters%n)), (/ 2, 2 /))

       !
       ! Include terms for speed
       !

       del2 = - del2 * hbar**2 / 2.0d0 / parameters%mass

       defdel = .true.

    end if

    if (defined_parameters%potential .and. (.not. defpot) .and.&
         & (.not. defvecpot)) then

       !
       ! Define the potential
       !

       allocate(pot(nl(1):nu(1), nl(2):nu(2)))

       call s_setpotential(pot, nl, nu)

       !
       ! Save the potential to disk for plotting
       !

       if (parameters%plotpot) then

          call s_write2drarray(real(pot), nl, nu, filename='rpot',&
               & incboundary=.true.)
          call s_write2drarray(aimag(pot), nl, nu, filename='ipot',&
               & incboundary=.true.)

       end if

       defpot = .true.

    end if

    if (defined_parameters%externalb .and. (.not. defvecpot) ) then

       select case (parameters%gauge)

       case('(-y,0,0)')

          allocate(delx(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfny(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn2(nl(1):nu(1), nl(2):nu(2)),&
               & delwavefnx(nl(1):nu(1), nl(2):nu(2)))

          call s_bpoty(delx, vecpotfny, vecpotfn2, nl, nu)

       case('(-y,x,0)')

          allocate(delx(nl(1):nu(1), nl(2):nu(2)), dely(nl(1):nu(1),&
               & nl(2):nu(2)), vecpotfnx(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfny(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn2(nl(1):nu(1), nl(2):nu(2)),&
               & delwavefnx(nl(1):nu(1), nl(2):nu(2)),&
               & delwavefny(nl(1):nu(1), nl(2):nu(2)))

          call s_bpotyx(delx, dely, vecpotfnx, vecpotfny, vecpotfn2,&
               & nl, nu)

       case('(x-y,x-y,0)')
          allocate(delxy(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn2(nl(1):nu(1), nl(2):nu(2)),&
               & delwavefn(nl(1):nu(1), nl(2):nu(2)))

          call s_bpotxyxy(delxy, vecpotfn, vecpotfn2, nl, nu)

       case('fft')
          allocate(delxy(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn(nl(1):nu(1), nl(2):nu(2)),&
               & vecpotfn2(nl(1):nu(1), nl(2):nu(2)),&
               & delwavefn(nl(1):nu(1), nl(2):nu(2)))

          call s_bpotfft(delxy, vecpotfn, vecpotfn2, nl, nu)

       end select

       defvecpot = .true.

       if (defpot) then

          vecpotfn2 = vecpotfn2 + pot

          deallocate(pot)

          defpot = .false.

       end if

    end if

    !
    ! Compute the derivatives of the wavefn
    !

    hamiltonian = wavefn

    call s_c2dfft(hamiltonian, nl, nu, fftfor)

    if (defined_parameters%externalb) then

       select case (parameters%gauge)

       case('(-y,0,0)')
          delwavefnx = hamiltonian*delx

          call s_c2dfft(delwavefnx, nl, nu, fftbac)

       case('(-y,x,0)')
          delwavefnx = hamiltonian*delx
          delwavefny = hamiltonian*dely

          call s_c2dfft(delwavefnx, nl, nu, fftbac)
          call s_c2dfft(delwavefny, nl, nu, fftbac)

       case('(x-y,x-y,0)', 'fft')
          delwavefn = hamiltonian*delxy

          call s_c2dfft(delwavefn, nl, nu, fftbac)

       end select

    end if

    hamiltonian = hamiltonian*del2

    call s_c2dfft(hamiltonian, nl, nu, fftbac)

    !
    ! Build up the hamiltonian from its various bits
    !

    if (defined_parameters%externalb .and. defvecpot) then

       select case (parameters%gauge)

       case('(-y,0,0)')
          hamiltonian = hamiltonian + vecpotfny*delwavefnx +&
               & vecpotfn2*wavefn

       case('(-y,x,0)')
          hamiltonian = hamiltonian + vecpotfnx*delwavefny +&
               & vecpotfny*delwavefnx + vecpotfn2*wavefn

       case('(x-y,x-y,0)', 'fft')
          hamiltonian = hamiltonian + vecpotfn*delwavefn + vecpotfn2&
               &*wavefn

       end select

    end if

    if (defined_parameters%potential .and. defpot) then

       hamiltonian = hamiltonian + pot*wavefn

    end if

    !
    ! If we have a magnetic field, ensure the boundary of the wave
    ! function
    ! is zero.  This is a nasty hack, but I think required.
    ! Hopefully it
    ! will make this work well for magnetic fields.
    !

    if(defined_parameters%externalb) then

       call s_smoothedge(hamiltonian, nl, nu, parameters%smooth)

    end if

  end function f_hamiltonian

  function f_wavefnenergy(wavefn, nl, nu) result(energy)

    !
    ! this module calculates the energy of a function
    !

    integer(idp), dimension(dim), intent(in) :: nl, nu
    complex(cdp), dimension(nl(1):nu(1), nl(2):nu(2)), intent(in) ::&
         & wavefn
    real(fdp) :: energy

    !
    !  Compute the hamiltonian
    !

    energy = f_norm(real(conjg(wavefn)*f_hamiltonian(wavefn, nl, nu))&
         &, nl, nu)

  end function f_wavefnenergy

end module hamiltonian