!> Module to handle a harmonic-oscillator model for the
!! bath in quantum dynamics simulations
MODULE harmonic_bath
  USE kinds
  IMPLICIT NONE
  REAL(dp), ALLOCATABLE :: x_bath(:), p_bath(:) ! bath phase space DOFs
  REAL(dp), ALLOCATABLE :: omega(:), c(:,:) ! bath frequencies and coupling matrix elements
  REAL(dp), ALLOCATABLE :: site_sd(:,:), site_reorg(:)
  REAL(dp), ALLOCATABLE :: site_bath_freq(:)
CONTAINS

    !> Allocates the necessary arrays to describe the bath and system-bath interactions. Also, the
    !! spectral density is uniformly sampled according to frequency-dependent reorganization energy.
  SUBROUTINE initialize_bath
    USE kinds
    USE input_output
    USE parameters
    USE unit_conversions
    IMPLICIT NONE
    INTEGER :: i,j,k,l,imode
    INTEGER :: npts ! Number of points in spectral density file
    INTEGER :: smode ! Sampled mode index
    CHARACTER :: junk ! Dummy character used to navigate input files
    REAL(dp) :: dw ! Discrete frequency element i.e.) freq(2) - freq(1)
    REAL(dp), ALLOCATABLE :: freq(:) ! Frequencies from continuumsd.in file
    REAL(dp), ALLOCATABLE :: sd(:,:), lambda(:,:) ! Spectral density from continuumsd.in, reorg engy
    REAL(dp), ALLOCATABLE :: norm(:) ! Normalization factor for each bath
    REAL(dp), ALLOCATABLE :: cfactor(:,:) ! Factor to turn on/off site<->bath coupling
    REAL(dp), ALLOCATABLE :: sampledw(:,:), sampledc(:,:), sampledj(:,:) ! Samplings from continuumsd.in
    dw = 0.0_dp
    !> Note: This is restricted to a single range and discretization
    !!       of frequencies (defined by a single dw)

    ! Figure out how many lines in spectral density file
    npts = 0
    OPEN(UNIT=10, FILE=SPECTRALIN)
    DO
        READ(10, *, ERR=11, END=11) junk
        npts=npts+1
    END DO
    11  CONTINUE ; CLOSE(10)

    ! Warn if there was no continuum spectral density provided
    IF (npts==0) THEN
        OPEN(unit=20,file=ERRORLOG)
        WRITE(20,*) 'No continuum spectral density file found!'
        CLOSE(20)
    END IF

    ! Figure out how many baths and how many oscillators per bath
    OPEN(UNIT=30,FILE=COUPLINGS)
    READ(30,*) nbath, junk, nosc
    READ(30,*) junk

    ! Allocate all arrays needed to describe the bath
    ALLOCATE(freq(npts), sd(nbath, npts), lambda(nbath, npts), norm(nbath), &
            cfactor(nbath, nstate), sampledw(nbath, nosc), sampledc(nbath, nosc), &
            sampledj(nbath, nosc), site_sd(nbath, npts), &
            site_reorg(nbath), site_bath_freq(npts))

    freq = 0.0_dp ; sd = 0.0_dp ; lambda = 0.0_dp ; cfactor = 0.0_dp
    sampledw = 0.0_dp ; sampledc = 0.0_dp ; sampledj = 0.0_dp ; norm = 0.0_dp
    site_sd = 0.0_dp ; site_reorg = 0.0_dp
    site_bath_freq = 0.0_dp

    ! Read in the prefactor matrix to turn on/off site<->bath coupling
    DO i=1, nstate
       READ(30,*) ( cfactor(j,i), j=1, nbath )
    END DO
    CLOSE(30)

    ! Read in the frequency and spectral density for each bath
    OPEN(UNIT=10,FILE=SPECTRALIN)
    DO i=1, npts
        READ(10, *) freq(i), ( sd(j,i), j=1, nbath )
    END DO
    CLOSE(10)

    ! Convert the frequencies and spectral densities to atomic units
    freq = freq * convert('wvnbr','au_ang_freq')
    sd = sd * convert('wvnbr','au_energy')

    dw = freq(2) - freq(1)

    site_bath_freq = freq
    site_sd = sd

    !> Sample bath according to the method described by
    !! Wang et al JCP 110, 4828 (1999)
    IF ( freq(1) /= 0.0_dp ) lambda(:,1) = dw * sd(:,1) / freq(1)
    DO i=1, nbath
       DO j=2, npts
          lambda(i,j) = lambda(i,j-1) + dw * sd(i,j) / freq(j)
       END DO
       site_reorg(i) = lambda(i,npts)
       norm(i) = DBLE(nosc) / lambda(i,npts)
       lambda(i,:) = lambda(i,:) * norm(i)
    END DO

   ! Start Sampling
    DO i=1, nbath
       smode = 1
       DO imode=1, nosc-1
          DO WHILE ( lambda(i,smode) <= REAL(imode) )
             smode = smode + 1
          END DO
          sampledw(i,imode) = freq(smode)
          sampledc(i,imode) = DSQRT( 2.0_dp / pi / norm(i) ) * freq(smode)
          sampledj(i,imode) = sd(i,smode)
       END DO
       sampledw(i,nosc) = freq(npts)
       sampledc(i,nosc) = DSQRT( 2.0_dp / pi / norm(i) ) * freq(npts)
       sampledj(i,nosc) = sd(i,npts)
    END DO

    ! Allocate coupling matrix, frequency vector, Bath DOF vectors
    ALLOCATE(x_bath(nbath*nosc), p_bath(nbath*nosc), omega(nbath*nosc), c(nbath*nosc, nstate))
    c = 0.0_dp ; omega = 0.0_dp ; x_bath = 0.0_dp ; p_bath = 0.0_dp

    ! Fill coupling and frequency arrays
    DO i=1, nstate
      DO l=1, nbath
         DO j=1, nosc
            k = ( l - 1 ) * nosc + j
            c(k,i) = cfactor(l,i) * sampledc(l,j)
            omega(k) = sampledw(l,j)
         END DO
      END DO
    END DO

   !> Writes a log of sampled modes into a file called "sampledbath.log"
   OPEN(unit=40,FILE=SPECTRALOUT)
   DO j=1, nosc
      WRITE(40,*) 'Freq', ( sampledw(i,j) * convert('au_ang_freq','wvnbr'), i = 1, nbath ) , &
           'Spec', ( sampledj(i,j) * convert('au_energy','wvnbr'), i=1, nbath )
   END DO
   CLOSE(40)

   DEALLOCATE(freq, sd, lambda, norm, cfactor, sampledw, sampledc, sampledj)

 END SUBROUTINE initialize_bath

    !> Sample the initial thermal distribution for harmonic bath from the
    !! Wigner transform of the Boltzmann operator. This is under the assumption that
    !! the full initil density operator can be written as \f$ \hat{\rho}_s\hat{\rho}_b^{eq} \f$.
 SUBROUTINE sample_thermal_wigner(x_init, p_init, beta)
   USE random_numbers
   USE unit_conversions
   USE input_output
   USE kinds
   IMPLICIT NONE
   REAL(dp), INTENT(in) :: beta ! For Boltzmann operator
   REAL(dp), INTENT(out) :: x_init(nbath*nosc), p_init(nbath*nosc) ! Initial bath DOFs

   ! Calcuate width of distribution for x_init, p_init
   p_init(:) = DSQRT( omega(:) / ( 2.0_dp * DTANH( 0.5_dp * beta * omega(:) ) ) )
   x_init(:) = p_init(:) / omega(:)

   ! Randomly sample the Distribution for x_init, p_init
   p_init(:) = p_init(:) * gaussian_rn(p_init)
   x_init(:) = x_init(:) * gaussian_rn(x_init)

 END SUBROUTINE sample_thermal_wigner

 ! Deallocate the bath arrays
 SUBROUTINE finalize_bath
   IMPLICIT NONE

   DEALLOCATE(x_bath, p_bath, omega, c)

 END SUBROUTINE finalize_bath

 !> Calculate bilinear coupling force for PLDM
 FUNCTION bilinear_harmonic_force_pldm(x, x_map, p_map, xt_map, pt_map) RESULT(res)
   USE kinds
   USE input_output
   IMPLICIT NONE
   INTEGER :: i
   REAL(dp) :: res(nbath*nosc)
   REAL(dp), INTENT(in) :: x(nbath*nosc)
   REAL(dp), INTENT(in) :: x_map(nstate), p_map(nstate), xt_map(nstate), pt_map(nstate)
   res(:) = - omega(:) ** 2 * x(:)
   DO i=1, nstate
      res(:) = res(:) - c(:,i) * 0.25_dp * ( x_map(i)**2 + p_map(i)**2 + xt_map(i)**2 + pt_map(i)**2 )
      ! 0.5 from the mapping variables and 0.5 from average of fwd & bkwd
   END DO

 END FUNCTION bilinear_harmonic_force_pldm

    !> Calculate the bilinear coupling force for TWA and SQC calculations
 FUNCTION bilinear_harmonic_force_twa(x, x_map, p_map) RESULT(res)
   USE kinds
   USE input_output
   IMPLICIT NONE
   INTEGER :: i
   REAL(dp) :: res(nbath*nosc)
   REAL(dp), INTENT(in) :: x(nbath*nosc)
   REAL(dp), INTENT(in) :: x_map(nstate), p_map(nstate)

   res(:) = - omega(:) ** 2 * x(:)
   DO i=1, nstate
      res(:) = res(:) - c(:,i) * 0.5_dp * ( x_map(i)**2 + p_map(i)**2 - 2.0_dp * zpe )
      ! 0.5 from the mapping variables
   END DO

 END FUNCTION bilinear_harmonic_force_twa

    !> Calculate the harmonic bath correlation function for a bilinearly coupled
    !! system-bath model in the Redfield Approximation for 1 ps
    FUNCTION redfield_bath_correlation_function(sd,w,beta_local) RESULT(res)
        USE kinds
        USE parameters
        USE input_output
        IMPLICIT NONE
        INTEGER :: itime, ifreq
        REAL(dp), INTENT(in) :: sd(:), w(:)
        REAL(dp), INTENT(in) :: beta_local
        REAL(dp) :: res(0 : nbstep), t, dw, n_of_w
        res = 0.0_dp ; dw = w(2) - w(1) ; n_of_w = 0.0_dp

        DO ifreq = 1, SIZE(w)
            n_of_w = 1.0_dp / ( exp(beta_local * w(ifreq)) - 1.0_dp )
            DO itime = 0, nbstep
                t = itime * ( runtime / DBLE(nbstep) )

                res(itime) = res(itime) + sd(ifreq) &
                        * ( ( n_of_w + 1.0_dp )* exp(-eye * w(ifreq) * t ) &
                + n_of_w * exp( eye * w(ifreq) * t ) )

            END DO
        END DO
        res = res * dw / pi
    END FUNCTION redfield_bath_correlation_function

    !> Subroutine to prepare the spectral density that is coupled
    !! to an "intramolecular" vibrational mode
    SUBROUTINE bath_coupled_to_vibrational_mode( q_init, pq_init, beta)
        USE input_output
        USE parameters
        USE kinds
        USE linear_algebra
        USE unit_conversions
        IMPLICIT NONE
        REAL(dp), INTENT(inout) :: q_init(:), pq_init(:)
        REAL(dp), INTENT(in) :: beta

    END SUBROUTINE bath_coupled_to_vibrational_mode

END MODULE harmonic_bath


!+---------+---------+---------+---------+---------+---------+--------+!
! this subroutine is used to calculate surface state using             !
! green's function method  ---  J.Phys.F.Met.Phys.15(1985)851-858      !
!
! History:
!         by Quan Sheng Wu on 4/20/2010                                !
!            mpi version      4/21/2010
!         Ca3PbO version     11/12/2011
!+---------+---------+---------+---------+---------+---------+--------+!
  subroutine fermiarc

     use wmpi
     use para
     implicit none

     integer :: ierr

     integer :: arclfile, arcrfile, arcbulkfile

     ! general loop index
     integer :: i,j,io

     ! kpoint loop index
     integer :: ikp

     real(dp) :: k(2)

     real(dp) :: omega
     real(dp) :: k1min_shape, k1max_shape, k2min_shape, k2max_shape

     real(dp), allocatable :: k12(:,:)
     real(dp), allocatable :: k12_shape(:,:)

     real(dp), allocatable :: dos_l(:)
     real(dp), allocatable :: dos_l_mpi(:)
     real(dp), allocatable :: dos_r(:)
     real(dp), allocatable :: dos_r_mpi(:)
     real(dp), allocatable :: dos_l_only(:)
     real(dp), allocatable :: dos_r_only(:)
     real(dp), allocatable :: dos_bulk(:)
     real(dp), allocatable :: dos_bulk_mpi(:)

     complex(dp), allocatable :: ones(:,:)
     complex(dp), allocatable :: GLL(:,:), GRR(:,:), GB(:,:)
     complex(dp), allocatable :: H00(:,:), H01(:,:)

     allocate( k12(2, nk1*nk2))
     allocate( k12_shape(2, nk1*nk2))
     allocate( dos_l(nk1*nk2))
     allocate( dos_l_mpi(nk1*nk2))
     allocate( dos_r(nk1*nk2))
     allocate( dos_l_only(nk1*nk2))
     allocate( dos_r_only(nk1*nk2))
     allocate( dos_r_mpi(nk1*nk2))
     allocate( dos_bulk(nk1*nk2))
     allocate( dos_bulk_mpi(nk1*nk2))
     allocate( GLL(ndim,ndim), GRR(ndim,ndim), GB(ndim, ndim))
     k12=0d0
     k12_shape=0d0
     dos_l=0d0
     dos_l_mpi=0d0
     dos_r=0d0
     dos_r_mpi=0d0
     dos_bulk=0d0
     dos_bulk_mpi=0d0

     ikp=0
     do i= 1, nk1
        do j= 1, nk2
           ikp=ikp+1
           k12(:, ikp)=K2D_start+ (i-1)*K2D_vec1/dble(nk1-1) &
                      + (j-1)*K2D_vec2/dble(nk2-1)
           k12_shape(:, ikp)= k12(1, ikp)* Ka2+ k12(2, ikp)* Kb2
        enddo
     enddo

     k1min_shape= minval(k12_shape(1, :))
     k2min_shape= minval(k12_shape(2, :))
     k1max_shape= maxval(k12_shape(1, :))
     k2max_shape= maxval(k12_shape(2, :))

     allocate(H00(Ndim, Ndim))
     allocate(H01(Ndim, Ndim))
     allocate(ones(Ndim, Ndim))
     GLL= 0d0
     GRR= 0d0
     H00= 0d0
     H01= 0d0
     ones= 0d0

     do i=1,Ndim
        ones(i,i)=1.0d0
     enddo

     omega = E_arc
     eta= eta_arc

     do ikp= 1+cpuid, nk1*nk2, num_cpu
        if (cpuid==0) write(stdout, *) 'Arc, ik, knv2', ikp, nk1*nk2
        k(1)= k12(1, ikp)
        k(2)= k12(2, ikp)

        call ham_qlayer2qlayer(k,H00,H01)

        !> calculate surface green function
        ! there are two method to calculate surface green's function
        ! the method in 1985 is better, you can find the ref in the
        ! subroutine
        call surfgreen_1985(omega,GLL,GRR,GB,H00,H01,ones)
        ! call surfgreen_1984(omega,GLL,GRR,H00,H01,ones)

        ! calculate spectral function
        do i= 1, NtopOrbitals
           io= TopOrbitals(i)
           dos_l(ikp)=dos_l(ikp)- aimag(GLL(io,io))
        enddo ! i
        do i= 1, NBottomOrbitals
           io= Ndim- Num_wann+ BottomOrbitals(i)
           dos_r(ikp)=dos_r(ikp)- aimag(GRR(io,io))
        enddo ! i
        do i= 1, Ndim
           dos_bulk(ikp)=dos_bulk(ikp)- aimag(GB (i,i))
        enddo ! i

     enddo

     !> we don't have to do allreduce operation
#if defined (MPI)
     call mpi_reduce(dos_l, dos_l_mpi, size(dos_l),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
     call mpi_reduce(dos_r, dos_r_mpi, size(dos_r),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
     call mpi_reduce(dos_bulk, dos_bulk_mpi, size(dos_bulk), mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
#else
     dos_l_mpi= dos_l
     dos_r_mpi= dos_r
     dos_bulk_mpi= dos_bulk
#endif

     do ikp=1, Nk1*Nk2
        dos_l_only(ikp)= dos_l_mpi(ikp)- dos_bulk_mpi(ikp)
        if (dos_l_only(ikp)<0) dos_l_only(ikp)=eps9
        dos_r_only(ikp)= dos_r_mpi(ikp)- dos_bulk_mpi(ikp)
        if (dos_r_only(ikp)<0) dos_r_only(ikp)=eps9
     enddo

     outfileindex= outfileindex+ 1
     arclfile= outfileindex
     outfileindex= outfileindex+ 1
     arcrfile= outfileindex
     outfileindex= outfileindex+ 1
     arcbulkfile= outfileindex
     if (cpuid.eq.0)then
        open (unit=arclfile, file='arc.dat_l')
        open (unit=arcrfile, file='arc.dat_r')
        open (unit=arcbulkfile, file='arc.dat_bulk')
        do ikp=1, nk1*nk2
           write(arclfile, '(30f16.8)')k12_shape(:, ikp), log(dos_l_mpi(ikp)), log(dos_l_only(ikp))
           if (mod(ikp, nk2)==0) write(arclfile, *)' '
           write(arcrfile, '(30f16.8)')k12_shape(:, ikp), log(dos_r_mpi(ikp)), log(dos_r_only(ikp))
           if (mod(ikp, nk2)==0) write(arcrfile, *)' '
           write(arcbulkfile, '(30f16.8)')k12_shape(:, ikp), log(dos_bulk_mpi(ikp))
           if (mod(ikp, nk2)==0) write(arcbulkfile, *)' '
        enddo
        close(arclfile)
        close(arcrfile)
        close(arcbulkfile)

        write(stdout,*)'ndim',ndim
        write(stdout,*)'Nk1,Nk2,eta', Nk1, Nk2, eta
        write(stdout,*)'calculate density of state successfully'
     endif

     !> write script for gnuplot for bulk green's function
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_bulk.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_bulk.png'"
        write(outfileindex,'(2a)') 'set palette defined ( -10 "#194eff", ', &
           '0 "white", 10 "red" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_bulk' u 1:2:3 w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_l_only.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_l_only.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_l' u 1:2:(exp($4)) w pm3d"

        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_l.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_l.png'"
        write(outfileindex,'(2a)') 'set palette defined ( -10 "#194eff", ', &
           '0 "white", 10 "red" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_l' u 1:2:3 w pm3d"

        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_r_only.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_r.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_r_only.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_r' u 1:2:(exp($4)) w pm3d"

        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_r.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_r.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_r.png'"
        write(outfileindex,'(2a)') 'set palette defined ( -10 "#194eff", ', &
           '0 "white", 10 "red" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_r' u 1:2:3 w pm3d"

        close(outfileindex)
     endif

  return
  end subroutine fermiarc

  !> calculation joint density of state
  !> refs :http://science.sciencemag.org/content/sci/suppl/2016/03/09/351.6278.1184.DC1/Inoue.SM.pdf

  !> in this version, we also calculate the spin density jdos
  subroutine fermiarc_jdos

     use wmpi
     use para
     implicit none

     integer :: ierr

     integer :: arclfile, arcrfile, arcbulkfile
     integer :: arcljfile, arcrjfile
     integer :: arcljsfile, arcrjsfile
     integer :: spindosrfile, spindoslfile

     ! general loop index
     integer :: i,j,io

     integer :: Nwann
     ! kpoint loop index
     integer :: ikp, ik1, ik2, iq, Nk1_half, Nk2_half
     integer :: imin1, imax1, imin2, imax2, iq1, iq2

     real(dp) :: k(2), dis

     real(dp) :: omega
     real(dp) :: k1min_shape, k1max_shape, k2min_shape, k2max_shape

     real(dp) :: sx_bulk, sy_bulk, sz_bulk

     integer , allocatable :: ik12(:,:)
     real(dp), allocatable :: k12(:,:)
     real(dp), allocatable :: k12_shape(:,:)

     real(dp), allocatable :: dos_l(:,:)
     real(dp), allocatable :: dos_l_mpi(:,:)
     real(dp), allocatable :: dos_r(:,:)
     real(dp), allocatable :: dos_r_mpi(:,:)
     real(dp), allocatable :: dos_bulk(:,:)
     real(dp), allocatable :: dos_bulk_mpi(:,:)
     real(dp), allocatable :: jdos_l(:)
     real(dp), allocatable :: jdos_l_mpi(:)
     real(dp), allocatable :: jdos_r(:)
     real(dp), allocatable :: jdos_r_mpi(:)
     real(dp), allocatable :: jsdos_l(:)
     real(dp), allocatable :: jsdos_l_mpi(:)
     real(dp), allocatable :: jsdos_r(:)
     real(dp), allocatable :: jsdos_r_mpi(:)
     real(dp), allocatable :: dos_l_only(:, :)
     real(dp), allocatable :: dos_r_only(:, :)
     real(dp), allocatable :: sx_l(:, :)
     real(dp), allocatable :: sy_l(:, :)
     real(dp), allocatable :: sz_l(:, :)
     real(dp), allocatable :: sx_l_mpi(:, :)
     real(dp), allocatable :: sy_l_mpi(:, :)
     real(dp), allocatable :: sz_l_mpi(:, :)
     real(dp), allocatable :: sx_r(:, :)
     real(dp), allocatable :: sy_r(:, :)
     real(dp), allocatable :: sz_r(:, :)
     real(dp), allocatable :: sx_r_mpi(:, :)
     real(dp), allocatable :: sy_r_mpi(:, :)
     real(dp), allocatable :: sz_r_mpi(:, :)

     ! spin operator matrix sigma_x,sigma_y in sigma_z representation
     complex(Dp),allocatable :: sigma_x(:,:)
     complex(Dp),allocatable :: sigma_y(:,:)
     complex(Dp),allocatable :: sigma_z(:,:)
     complex(Dp),allocatable :: ctemp(:,:)
     complex(dp), allocatable :: ones(:,:)
     complex(dp), allocatable :: GLL(:,:), GRR(:,:), GB(:,:)
     complex(dp), allocatable :: H00(:,:), H01(:,:)

     Nk1_half= (Nk1-1)/2
     Nk2_half= (Nk2-1)/2
     Nwann= Num_wann/2

     allocate( ik12(2, nk1*nk2))
     allocate( k12(2, nk1*nk2))
     allocate( k12_shape(2, nk1*nk2))
     allocate( dos_l(nk1,nk2))
     allocate( dos_l_mpi(nk1,nk2))
     allocate( dos_r(nk1,nk2))
     allocate( dos_r_mpi(nk1,nk2))
     allocate( dos_r_only(nk1,nk2))
     allocate( dos_l_only(nk1,nk2))
     allocate( dos_bulk(nk1,nk2))
     allocate( dos_bulk_mpi(nk1,nk2))
     allocate( jdos_l(nk1*nk2))
     allocate( jdos_l_mpi(nk1*nk2))
     allocate( jdos_r(nk1*nk2))
     allocate( jdos_r_mpi(nk1*nk2))
     allocate( jsdos_l(nk1*nk2))
     allocate( jsdos_l_mpi(nk1*nk2))
     allocate( jsdos_r(nk1*nk2))
     allocate( jsdos_r_mpi(nk1*nk2))
     allocate( sx_l(nk1,nk2))
     allocate( sx_l_mpi(nk1,nk2))
     allocate( sy_l(nk1,nk2))
     allocate( sy_l_mpi(nk1,nk2))
     allocate( sz_l(nk1,nk2))
     allocate( sz_l_mpi(nk1,nk2))
     allocate( sx_r(nk1,nk2))
     allocate( sx_r_mpi(nk1,nk2))
     allocate( sy_r(nk1,nk2))
     allocate( sy_r_mpi(nk1,nk2))
     allocate( sz_r(nk1,nk2))
     allocate( sz_r_mpi(nk1,nk2))
     allocate( GLL(ndim,ndim), GRR(ndim,ndim), GB(ndim, ndim))
     allocate(sigma_x(ndim,ndim))
     allocate(sigma_y(ndim,ndim))
     allocate(sigma_z(ndim,ndim))
     allocate(ctemp(ndim,ndim))

     sigma_x   =0.0d0
     sigma_y   =0.0d0
     sigma_z   =0.0d0
     ik12=0
     k12=0d0
     k12_shape=0d0
     dos_l=0d0
     dos_l_mpi=0d0
     dos_r=0d0
     dos_r_mpi=0d0
     jdos_l=0d0
     jdos_l_mpi=1d-12
     jdos_r=0d0
     jdos_r_mpi=1d-12
     jsdos_l=0d0
     jsdos_l_mpi=1d-12
     jsdos_r=0d0
     jsdos_r_mpi=1d-12
     sx_l=0d0
     sy_l=0d0
     sz_l=0d0
     sx_l_mpi=0d0
     sy_l_mpi=0d0
     sz_l_mpi=0d0
     sx_r=0d0
     sy_r=0d0
     sz_r=0d0
     sx_r_mpi=0d0
     sy_r_mpi=0d0
     sz_r_mpi=0d0

     ikp=0
     do i= 1, nk1
        do j= 1, nk2
           ikp=ikp+1
           ik12(1, ikp)= i
           ik12(2, ikp)= j
           k12(:, ikp)=K2D_start+ (i-1)*K2D_vec1/dble(nk1-1) &
                      + (j-1)*K2D_vec2/dble(nk2-1)
           k12_shape(:, ikp)= k12(1, ikp)* Ka2+ k12(2, ikp)* Kb2
        enddo
     enddo

     k1min_shape= minval(k12_shape(1, :))
     k2min_shape= minval(k12_shape(2, :))
     k1max_shape= maxval(k12_shape(1, :))
     k2max_shape= maxval(k12_shape(2, :))

     allocate(H00(Ndim, Ndim))
     allocate(H01(Ndim, Ndim))
     allocate(ones(Ndim, Ndim))
     GLL= 0d0
     GRR= 0d0
     H00= 0d0
     H01= 0d0
     ones= 0d0

     do i=1,Ndim
        ones(i,i)=1.0d0
     enddo

     Nwann= Num_wann/2

     !> spin operator matrix
     do i=1, Np
        do j=1, Nwann
           sigma_x(Num_wann*(i-1)+j, Num_wann*(i-1)+Nwann+j)=1.0d0
           sigma_x(Num_wann*(i-1)+j+Nwann, Num_wann*(i-1)+j)=1.0d0
           sigma_y(Num_wann*(i-1)+j, Num_wann*(i-1)+Nwann+j)=-zi
           sigma_y(Num_wann*(i-1)+j+Nwann, Num_wann*(i-1)+j)=zi
           sigma_z(Num_wann*(i-1)+j, Num_wann*(i-1)+j)= 1d0
           sigma_z(Num_wann*(i-1)+j+Nwann, Num_wann*(i-1)+j+Nwann)=-1d0
        enddo
     enddo

     omega = E_arc
     eta= eta_arc

     do ikp= 1+cpuid, nk1*nk2, num_cpu
        if (cpuid==0) write(stdout, *) 'Arc, ik, knv2', ikp, nk1*nk2
        k(1)= k12(1, ikp)
        k(2)= k12(2, ikp)

        call ham_qlayer2qlayer(k,H00,H01)

        !> calculate surface green function
        ! there are two method to calculate surface green's function
        ! the method in 1985 is better, you can find the ref in the
        ! subroutine
        call surfgreen_1985(omega,GLL,GRR,GB,H00,H01,ones)
        ! call surfgreen_1984(omega,GLL,GRR,H00,H01,ones)

        ik1= ik12(1, ikp)
        ik2= ik12(2, ikp)

        ! calculate spectral function
        do i= 1, NtopOrbitals
           io= TopOrbitals(i)
           dos_l(ik1, ik2)=dos_l(ik1, ik2)- aimag(GLL(io,io))
        enddo ! i
        do i= 1, NBottomOrbitals
           io= Ndim- Num_wann+ BottomOrbitals(i)
           dos_r(ik1, ik2)=dos_r(ik1, ik2)- aimag(GRR(io,io))
        enddo ! i

        do i= 1, Ndim
        !  dos_l(ik1, ik2)=dos_l(ik1, ik2)- aimag(GLL(i,i))
        enddo ! i
        do i= 1, Ndim
        !  dos_r(ik1, ik2)=dos_r(ik1, ik2)- aimag(GRR(i,i))
        enddo ! i
        do i= 1, Ndim
           dos_bulk(ik1, ik2)=dos_bulk(ik1, ik2)- aimag(GB (i,i))
        enddo ! i

        !>> calculate spin-resolved bulk spectrum
        sx_bulk= 0d0
        call mat_mul(ndim,GB ,sigma_x,ctemp)
        do i= 1, Ndim
           sx_bulk=sx_bulk- aimag(ctemp(i,i))
        enddo ! i

        sy_bulk= 0d0
        call mat_mul(ndim,GB ,sigma_y,ctemp)
        do i= 1, Ndim
           sy_bulk=sy_bulk- aimag(ctemp(i,i))
        enddo ! i

        sz_bulk= 0d0
        call mat_mul(ndim,GB ,sigma_z,ctemp)
        do i= 1, Ndim
           sz_bulk=sz_bulk- aimag(ctemp(i,i))
        enddo ! i

        !>> calculate spin-resolved surface spectrum

        !ctemp=matmul(surfgreen,sigma_x)
        call mat_mul(ndim,GLL,sigma_x,ctemp)
        do i= 1, NtopOrbitals
           io= TopOrbitals(i)
           sx_l(ik1, ik2)=sx_l(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sx_l(ik1, ik2)= sx_l(ik1, ik2)- sx_bulk
        if (sx_l(ik1, ik2)<0) sx_l(ik1, ik2)= eps9

        !ctemp=matmul(surfgreen,sigma_y)
        call mat_mul(ndim,GLL,sigma_y,ctemp)
        do i= 1, NtopOrbitals
           io= TopOrbitals(i)
           sy_l(ik1, ik2)=sy_l(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sy_l(ik1, ik2)= sy_l(ik1, ik2)- sy_bulk
        if (sy_l(ik1, ik2)<0) sy_l(ik1, ik2)= eps9

        !ctemp=matmul(surfgreen,sigma_z)
        call mat_mul(ndim,GLL,sigma_z,ctemp)
        do i= 1, NtopOrbitals
           io= TopOrbitals(i)
           sz_l(ik1, ik2)=sz_l(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sz_l(ik1, ik2)= sz_l(ik1, ik2)- sz_bulk
        if (sz_l(ik1, ik2)<0) sz_l(ik1, ik2)= eps9

        !ctemp=matmul(surfgreen,sigma_x)
        call mat_mul(ndim,GRR,sigma_x,ctemp)
        do i= 1, NBottomOrbitals
           io= Ndim- Num_wann+ BottomOrbitals(i)
           sx_r(ik1, ik2)=sx_r(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sx_r(ik1, ik2)= sx_r(ik1, ik2)- sz_bulk
        if (sx_r(ik1, ik2)<0) sx_r(ik1, ik2)= eps9

        !ctemp=matmul(surfgreen,sigma_y)
        call mat_mul(ndim,GRR,sigma_y,ctemp)
        do i= 1, NBottomOrbitals
           io= Ndim- Num_wann+ BottomOrbitals(i)
           sy_r(ik1, ik2)=sy_r(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sy_r(ik1, ik2)= sy_r(ik1, ik2)- sz_bulk
        if (sy_r(ik1, ik2)<0) sy_r(ik1, ik2)= eps9

        !ctemp=matmul(surfgreen,sigma_z)
        call mat_mul(ndim,GRR,sigma_z,ctemp)
        do i= 1, NBottomOrbitals
           io= Ndim- Num_wann+ BottomOrbitals(i)
           sz_r(ik1, ik2)=sz_r(ik1, ik2)- aimag(ctemp(io,io))
        enddo ! i
        sz_r(ik1, ik2)= sz_r(ik1, ik2)- sz_bulk
        if (sz_r(ik1, ik2)<0) sz_r(ik1, ik2)= eps9

     enddo

#if defined (MPI)
     !> we don't have to do allreduce operation
     call mpi_allreduce(dos_l, dos_l_mpi, size(dos_l),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(dos_r, dos_r_mpi, size(dos_r),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(dos_bulk, dos_bulk_mpi, size(dos_bulk),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sx_l, sx_l_mpi, size(sx_l),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sy_l, sy_l_mpi, size(sy_l),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sz_l, sz_l_mpi, size(sz_l),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sx_r, sx_r_mpi, size(sx_r),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sy_r, sy_r_mpi, size(sy_r),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)
     call mpi_allreduce(sz_r, sz_r_mpi, size(sz_r),mpi_double_precision,&
                     mpi_sum, mpi_comm_world, ierr)

#else
     dos_l_mpi= dos_l
     dos_r_mpi= dos_r
     dos_bulk_mpi= dos_bulk
     sx_l_mpi= sx_l
     sy_l_mpi= sy_l
     sz_l_mpi= sz_l
     sx_r_mpi= sx_r
     sy_r_mpi= sy_r
     sz_r_mpi= sz_r
#endif

     do ikp=1, Nk1*Nk2
        ik1= ik12(1, ikp)
        ik2= ik12(2, ikp)

        dos_l_only(ik1, ik2)= dos_l_mpi(ik1, ik2)- dos_bulk_mpi(ik1, ik2)
        if (dos_l_only(ik1, ik2)<0) dos_l_only(ik1, ik2)=eps9
        dos_r_only(ik1, ik2)= dos_r_mpi(ik1, ik2)- dos_bulk_mpi(ik1, ik2)
        if (dos_r_only(ik1, ik2)<0) dos_r_only(ik1, ik2)=eps9
     enddo

     outfileindex= outfileindex+ 1
     arclfile= outfileindex
     outfileindex= outfileindex+ 1
     arcrfile= outfileindex
     outfileindex= outfileindex+ 1
     spindoslfile= outfileindex
     outfileindex= outfileindex+ 1
     spindosrfile= outfileindex
     if (cpuid.eq.0)then
        open (unit=arclfile, file='arc.dat_l')
        open (unit=arcrfile, file='arc.dat_r')
        open(spindoslfile,file='spindos_l.dat')
        open(spindoslfile,file='spindos_r.dat')
        do ikp=1, nk1*nk2
           write(arclfile, '(30f16.8)')k12_shape(:, ikp), log(dos_l_mpi(ik12(1, ikp), ik12(2, ikp))) &
                                     , log(dos_l_only(ik12(1, ikp), ik12(2, ikp)))
           if (mod(ikp, nk2)==0) write(arclfile, *)' '
           write(arcrfile, '(30f16.8)')k12_shape(:, ikp), log(dos_r_mpi(ik12(1, ikp), ik12(2, ikp))) &
                                     , log(dos_r_only(ik12(1, ikp), ik12(2, ikp)))
           if (mod(ikp, nk2)==0) write(arcrfile, *)' '
           write(spindoslfile, '(30f16.8)')k12_shape(:, ikp), &
                                     log(dos_l_only(ik12(1, ikp), ik12(2, ikp))), &
                                     (sx_l_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_l_mpi(ik12(1, ikp), ik12(2, ikp)), &
                                     (sy_l_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_l_mpi(ik12(1, ikp), ik12(2, ikp)), &
                                     (sz_l_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_l_mpi(ik12(1, ikp), ik12(2, ikp))
           write(spindoslfile, '(30f16.8)')k12_shape(:, ikp), &
                                     log(dos_r_only(ik12(1, ikp), ik12(2, ikp))), &
                                     (sx_r_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_r_mpi(ik12(1, ikp), ik12(2, ikp)), &
                                     (sy_r_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_r_mpi(ik12(1, ikp), ik12(2, ikp)), &
                                     (sz_r_mpi(ik12(1, ikp), ik12(2, ikp)))/dos_r_mpi(ik12(1, ikp), ik12(2, ikp))
           if (mod(ikp, nk2)==0) write(spindoslfile, *)' '
           if (mod(ikp, nk2)==0) write(spindoslfile, *)' '
        enddo
        close(arclfile)
        close(arcrfile)
        close(spindoslfile)
        close(spindoslfile)
        write(stdout,*)'ndim',ndim
        write(stdout,*) 'Nk1,Nk2,eta',Nk1, Nk2, eta
        write(stdout,*)'calculate density of state successfully'
     endif

     !> calculate jdos
     do iq= 1+ cpuid, Nk1*Nk2, num_cpu
        iq1= ik12(1, iq)- Nk1_half
        iq2= ik12(2, iq)- Nk2_half
        if (cpuid==0) write(stdout, *) 'Jdos, iq, knv2', iq, nk1*nk2
        imin1= max(-Nk1_half-iq1, -Nk1_half)+ Nk1_half+ 1
        imax1= min(Nk1_half-iq1, Nk1_half)+ Nk1_half+ 1
        imin2= max(-Nk2_half-iq2, -Nk2_half)+ Nk2_half+ 1
        imax2= min(Nk2_half-iq2, Nk2_half)+ Nk2_half+ 1
        do ik1= imin1, imax1
           do ik2= imin2, imax2
              jdos_l(iq)= jdos_l(iq)+ dos_l_only(ik1, ik2)* dos_l_only(ik1+iq1, ik2+iq2)
              jdos_r(iq)= jdos_r(iq)+ dos_r_only(ik1, ik2)* dos_r_only(ik1+iq1, ik2+iq2)
              jsdos_l(iq)= jsdos_l(iq)+ dos_l_only(ik1, ik2)* dos_l_only(ik1+iq1, ik2+iq2) &
                                      + sx_l_mpi(ik1, ik2)* sx_l_mpi(ik1+iq1, ik2+iq2) &
                                      + sy_l_mpi(ik1, ik2)* sy_l_mpi(ik1+iq1, ik2+iq2) &
                                      + sz_l_mpi(ik1, ik2)* sz_l_mpi(ik1+iq1, ik2+iq2)
              jsdos_r(iq)= jsdos_r(iq)+ dos_r_only(ik1, ik2)* dos_r_only(ik1+iq1, ik2+iq2) &
                                      + sx_r_mpi(ik1, ik2)* sx_r_mpi(ik1+iq1, ik2+iq2) &
                                      + sy_r_mpi(ik1, ik2)* sy_r_mpi(ik1+iq1, ik2+iq2) &
                                      + sz_r_mpi(ik1, ik2)* sz_r_mpi(ik1+iq1, ik2+iq2)
           enddo !ik
        enddo !ik
     enddo !ik

     jdos_l_mpi=1d-12
     jdos_r_mpi=1d-12
     jsdos_l_mpi=1d-12
     jsdos_r_mpi=1d-12
#if defined (MPI)
     call mpi_reduce(jdos_l, jdos_l_mpi, size(jdos_l),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
     call mpi_reduce(jdos_r, jdos_r_mpi, size(jdos_r),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
     call mpi_reduce(jsdos_l, jsdos_l_mpi, size(jsdos_l),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)
     call mpi_reduce(jsdos_r, jsdos_r_mpi, size(jsdos_r),mpi_double_precision,&
                     mpi_sum, 0, mpi_comm_world, ierr)

#else
     jdos_l_mpi= jdos_l
     jdos_r_mpi= jdos_r
     jsdos_l_mpi= jsdos_l
     jsdos_r_mpi= jsdos_r
#endif

     outfileindex= outfileindex+ 1
     arcljfile= outfileindex
     outfileindex= outfileindex+ 1
     arcrjfile= outfileindex
     outfileindex= outfileindex+ 1
     arcljsfile= outfileindex
     outfileindex= outfileindex+ 1
     arcrjsfile= outfileindex
     if (cpuid.eq.0)then
        open (unit=arcljfile, file='arc.jdat_l')
        open (unit=arcrjfile, file='arc.jdat_r')
        open (unit=arcljsfile, file='arc.jsdat_l')
        open (unit=arcrjsfile, file='arc.jsdat_r')
        do ikp=1, nk1*nk2
           write(arcljfile, '(30f16.8)')k12_shape(:, ikp), log(jdos_l_mpi(ikp))
           if (mod(ikp, nk2)==0) write(arcljfile, *)' '
           write(arcrjfile, '(30f16.8)')k12_shape(:, ikp), log(jdos_r_mpi(ikp))
           if (mod(ikp, nk2)==0) write(arcrjfile, *)' '
           write(arcljsfile, '(30f16.8)')k12_shape(:, ikp), log(jsdos_l_mpi(ikp))
           if (mod(ikp, nk2)==0) write(arcljsfile, *)' '
           write(arcrjsfile, '(30f16.8)')k12_shape(:, ikp), log(jsdos_r_mpi(ikp))
           if (mod(ikp, nk2)==0) write(arcrjsfile, *)' '
        enddo
        close(arcljsfile)
        close(arcrjsfile)
        close(arcljsfile)
        close(arcrjsfile)
        write(stdout,*)'calculate joint density of state successfully'
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_l_jsdos.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l_jsdos.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_l_jsdos.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.jsdat_l' u 1:2:(exp($3)) w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_l_jdos.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l_jdos.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_l_jdos.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.jdat_l' u 1:2:(exp($3)) w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_l_only.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_l_only.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_l_only.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_l' u 1:2:(exp($4)) w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_r_jsdos.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_r_jsdos.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_r_jsdos.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.jsdat_r' u 1:2:(exp($3)) w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_r_jdos.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_r_jdos.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_r_jdos.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.jdat_r' u 1:2:(exp($3)) w pm3d"
        close(outfileindex)
     endif

     !> write script for gnuplot
     outfileindex= outfileindex+ 1
     if (cpuid==0) then
        open(unit=outfileindex, file='arc_r_only.gnu')
        write(outfileindex, '(a)')"set encoding iso_8859_1"
        write(outfileindex, '(a)')'#set terminal  postscript enhanced color'
        write(outfileindex, '(a)')"#set output 'arc_r_only.eps'"
        write(outfileindex, '(3a)')'#set terminal pngcairo truecolor enhanced', &
           '  font ",50" size 1920, 1680'
        write(outfileindex, '(3a)')'set terminal png truecolor enhanced', &
           ' font ",50" size 1920, 1680'
        write(outfileindex, '(a)')"set output 'arc_r_only.png'"
        write(outfileindex,'(2a)') 'set palette defined (0  "white", ', &
           '6 "red", 20 "black" )'
        write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10'
        write(outfileindex, '(a)')'unset ztics'
        write(outfileindex, '(a)')'unset key'
        write(outfileindex, '(a)')'set pm3d'
        write(outfileindex, '(a)')'set border lw 6'
        write(outfileindex, '(a)')'set size ratio -1'
        write(outfileindex, '(a)')'set view map'
        write(outfileindex, '(a)')'set xtics'
        write(outfileindex, '(a)')'set ytics'
        write(outfileindex, '(a)')'set xlabel "K_1 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel "K_2 (1/{\305})"'
        write(outfileindex, '(a)')'set ylabel offset 1, 0'
        write(outfileindex, '(a)')'set colorbox'
        write(outfileindex, '(a)')'unset cbtics'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set xrange [', k1min_shape, ':', k1max_shape, ']'
        write(outfileindex, '(a, f8.5, a, f8.5, a)')'set yrange [', k2min_shape, ':', k2max_shape, ']'
        write(outfileindex, '(a)')'set pm3d interpolate 2,2'
        write(outfileindex, '(2a)')"splot 'arc.dat_r' u 1:2:(exp($4)) w pm3d"

        close(outfileindex)
     endif

  return
  end subroutine fermiarc_jdos


text:
C$Procedure      SHIFTL ( Shift left )

      SUBROUTINE SHIFTL ( IN, NSHIFT, FILLC, OUT )

C$ Abstract
C
C      Shift the contents of a character string to the left.
C      Characters moved past the beginning of the input string are
C      lost. Vacant spaces are filled with a specified character.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C      CHARACTER,  UTILITY
C
C$ Declarations

      CHARACTER*(*)      IN
      INTEGER            NSHIFT
      CHARACTER*1        FILLC
      CHARACTER*(*)      OUT

C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C      IN         I   Input string.
C      NSHIFT     I   Number of times to shift.
C      FILLC      I   Character to fill spaces left vacant.
C      OUT        O   Shifted string.
C
C$ Detailed_Input
C
C      IN          is the input character string.
C
C      NSHIFT      is the number of times the string is to be
C                  shifted. If NSHIFT is negative, OUT will be
C                  identical to IN.
C
C      FILLC       is the character with which spaces left vacant by
C                  the shift are to be filled.
C
C$ Detailed_Output
C
C      OUT         is the output string. This is the input string,
C                  shifted N times, filled with FILLC.
C
C                  OUT may overwrite IN.
C
C$ Parameters
C
C     None.
C
C$ Particulars
C
C      As a string is shifted left or right, the leftmost or
C      rightmost characters of the string disappear (as if pushed
C      off the end of the string). This is true regardless of
C      the length of the output string.
C
C      The remaining characters are shifted simultaneously, and
C      the spaces vacated by those characters are filled with a
C      replacement character.
C
C$ Examples
C
C    If FILLC = ' '
C
C       'abcde'   shifted left twice becomes     'cde  '
C       'abcde'   shifted right once becomes     ' abcd'
C
C    If FILLC = '.'
C
C       '12345 '  shifted right once becomes     '.12345'
C       'Apple '  shifted left ten times becomes '......'
C
C    Given the declarations
C
C       CHARACTER*3         SHORT
C       CHARACTER*10        LONG
C
C    The calls
C
C       CALL SHIFTR ( 'abcde ', 2, '-', SHORT )
C       CALL SHIFTR ( 'abcde ', 2, '-', LONG  )
C
C    yield the strings
C
C       SHORT = '--a'
C       LONG  = '--abcd    '
C
C    while the calls
C
C       CALL SHIFTL ( 'abcde ', 2, '-', SHORT )
C       CALL SHIFTL ( 'abcde ', 2, '-', LONG  )
C
C    yield the strings
C
C       SHORT = 'cde'
C       LONG  = 'cde ..     '
C
C$ Restrictions
C
C      None.
C
C$ Exceptions
C
C      Error free.
C
C$ Files
C
C      None.
C
C$ Author_and_Institution
C
C      M.J. Spencer    (JPL)
C      I.M. Underwood  (JPL)
C
C$ Literature_References
C
C      None.
C
C$ Version
C
C-     SPICELIB Version 2.0.1, 22-AUG-2001 (EDW)
C
C         Corrected ENDDO to END DO.
C
C-     SPICELIB Version 2.0.0, 01-SEP-1994 (MJS)
C
C         This version correctly handles negative shifts.
C
C-     SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990 (IMU)
C
C-&

C$ Index_Entries
C
C     shift left
C
C-&

C
C     Local variables
C
      INTEGER          INLEN
      INTEGER          OUTLEN

      INTEGER          N
      INTEGER          NSAVE
      INTEGER          NFILL
      INTEGER          I
      INTEGER          S

C
C     Get the length of the input, output strings.
C
      INLEN  = LEN ( IN  )
      OUTLEN = LEN ( OUT )

C
C     If the shift is zero or negative, the string is not changed.
C     If longer than the input string, the entire string is shifted.
C
      S = MAX ( NSHIFT, 0 )
      N = MIN ( INLEN,  S )

C
C     Figure out how many characters in the input string will
C     be saved (will not be shifted off the end of the string,
C     and will fit in the output string), and how many fill
C     characters will be needed (no more than NSHIFT, no fewer
C     than zero).
C
      NSAVE =     MIN ( INLEN - N, OUTLEN         )
      NFILL = N - MAX ( 0,         INLEN - OUTLEN )

C
C     Move the saved characters to output.
C
      DO I = 1, NSAVE
         OUT(I:I) = IN(I+S:I+S)
      END DO

C
C     Add as many fill characters as appropriate.
C
      DO I = NSAVE+1, NSAVE+NFILL
         OUT(I:I) = FILLC
      END DO

C
C     Pad the output string with blanks (to cover any previous
C     ugliness there).
C
      IF ( OUTLEN .GT. INLEN ) THEN
         OUT(INLEN+1: ) = ' '
      END IF

      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


! This module stores the runtime parameters.
! These parameter are initialized in set_method_params().

module meth_params_module

  use amrex_fort_module, only : rt => amrex_real

  implicit none

  real(rt), save :: difmag        ! used only in consup to weight the divu contributin
  integer , save :: iorder        ! used only in uslope and uflaten

  real(rt), save, public  :: gamma_const, gamma_minus_1

  integer, parameter     :: NHYP    = 4
  integer, parameter     :: MAXADV  = 5

  ! NTHERM: number of thermodynamic variables
  integer         , save :: NTHERM, NVAR, NDIAG
  integer         , save :: URHO, UMX, UMY, UMZ, UEDEN, UEINT, UFA, UFS, UFX
  integer         , save :: TEMP_COMP, NE_COMP, ZHI_COMP

  ! QTHERM: number of primitive variables
  integer         , save :: QTHERM, QVAR
  integer         , save :: QRHO, QU, QV, QW, QPRES, QREINT, QFA, QFS

  integer         , save :: nadv

  real(rt)        , save :: small_dens, small_temp, small_pres

  integer         , save :: ppm_type
  integer         , save :: ppm_reference
  integer         , save :: ppm_flatten_before_integrals
  integer         , save :: use_colglaz
  integer         , save :: use_flattening
  integer         , save :: version_2
  integer         , save :: corner_coupling
  integer         , save :: use_const_species
  integer         , save :: normalize_species
  integer         , save :: heat_cool_type
  integer         , save :: inhomo_reion
  integer         , save :: grav_source_type

  ! indices for the axion field
  integer         , save :: NAXVAR
  integer         , save :: UAXDENS, UAXRE, UAXIM, UAXPHAS

  integer, save :: npassive
  integer, save, allocatable :: qpass_map(:), upass_map(:)

end module meth_params_module


text:
Program XrayTrans685

Use Kinds, Only: dp
Use Kinds, Only: id
Use FileIO_Utilities, Only: Delta_Time
Implicit None
Integer(id) :: c

Call SYSTEM_CLOCK(COUNT=c)
Call Transport_Photons( '1e_5__' , 1000000 , hv0_pdf_Line_5 , (/1.E1_dp,1.E2_dp,1.E3_dp,1.E4_dp,1.E5_dp/) )
Write(*,'(ES12.3E2,A)') Real(60_id*1000000_id,dp)/Delta_Time(clock_then=c),' histories per minute'

Contains

Subroutine Transport_Photons(p_name,n_p,hv0_pdf,Rd)
    Use Kinds, Only: dp
    Use Global, Only: default_seed
    Use Global, Only: speed_of_light_SI
    Use Global, Only: m0c2
    Use Global, Only: Z_hat,Y_hat
    Use Global, Only: pi
    Use FileIO_Utilities, Only: creturn
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Character, Intent(In) :: p_name  !problem name, will start output filenames
    Integer, Intent(In) :: n_p  !number of particles to be run
    Interface
        Function hv0_pdf(r) Result(hv0)
            Use Kinds, Only: dp
            Use PRNGs, Only: ACORN_Type
            Real(dp) :: hv0
            Type(ACORN_Type), Intent(InOut) :: r
        End Function hv0_pdf
    End Interface
    Real(dp), Intent(In) :: Rd(:)  ![cm] radii of detector shells

    Real(dp), Parameter :: rho0 = 1225.0002517923321962E-6_dp  ![g/cm^3] atmospheric density at sea level
    Real(dp), Parameter :: speed_c = speed_of_light_SI * 100._dp  ![cm/s] speed of light
    Integer, Parameter :: photoelectric_absorption = 1
    Integer, Parameter :: incoherent_scatter = 2
    Integer, Parameter :: coherent_scatter = 3
    Integer, Parameter :: pair_production = 4

    Type(ACORN_Type) :: RNG
    Integer :: n_sig  !number of lines in the MAC list
    Real(dp), Allocatable :: hv_sigs(:)  ![keV] energy list for the MAC list, stored as logarithms
    Real(dp), Allocatable :: sigs(:,:)  ![cm^2/g] MAC lists, stored as logarithms
    Integer, Parameter :: n_E = 80
    Real(dp) :: E_list(0:n_E)
    Integer :: n_d
    Integer :: p_list_size
    Real(dp), Allocatable :: p_list(:,:,:)
    Integer, Allocatable :: p_listed(:)
    Real(dp), Allocatable :: f(:,:)
    Real(dp), Allocatable :: f_list(:,:,:)
    Integer :: p  !counter for particles
    Real(dp) :: r(1:3)  ![cm] position of the current photon
    Real(dp) :: omega_hat(1:3)  ![1 cm/s] direction of travel of the current photon
    Real(dp) :: hv  ![keV] energy of the current photon
    Real(dp) :: t  ![s] total flight time of the current photon
    Real(dp) :: w  !weight of the current photon
    Real(dp) :: r0  ![cm] magnitude of the photon position vector, before the most recent movement
    Real(dp) :: zeta0  !cosine of the zenith angle for the direction of travel for the current photon
    Real(dp) :: tau  ![mfp] optical thickness to the next interaction
    Real(dp) :: sig(1:4)  ![cm^2/g] mass attenuation coefficients for absorption, incoherent scat, coherent scat, and pair prod
    Real(dp) :: sigT  ![cm^2/g] total mass attenuation coefficient
    Real(dp) :: d  ![cm] geometric distance to the next interaction
    Real(dp) :: r1
    Real(dp) :: zeta
    Real(dp) :: Rca
    Real(dp) :: dca
    Real(dp) :: dz
    Real(dp) :: xi  !a scaled random number
    Integer :: interaction  !index of the next interaction
    Real(dp) :: sig_sum  !running sum of mass attenuation coefficients
    Real(dp) :: alpha  !ratio of photon energy to rest mass energy of an electron
    Real(dp) :: mu0 !cosine of the polar scattering angle
    Real(dp) :: w0  ![rad] azimuthal scattering angle
    Real(dp) :: b_hat(1:3)  ![1 cm] basis vector for correlating scattered direction to incident direction
    Real(dp) :: c_hat(1:3)  ![1 cm] basis vector for correlating scattered direction to incident direction
    Real(dp) :: nu  !sine of the polar scattering angle
    Real(dp) :: fp
    Integer :: i,j,s

    !Initialize random number generator
    Call RNG%seed(seed=default_seed)
    !Read in MAC values
    Call Read_MAC('resources/mac/07-N.txt',n_sig,hv_sigs,sigs)
    !Set up detector energy grid
    Do i = 0,n_E
        E_list(i) = 10._dp**(-4._dp + 0.1_dp*Real(i,dp))
    End Do
    !Initialize tally arrays
    n_d = Size(Rd)
    p_list_size = 2*n_d
    Allocate(p_list(1:2,1:p_list_size,1:n_d))
    p_list = 0._dp
    Allocate(p_listed(1:n_d))
    p_listed = 0
    Allocate(f(1:2,1:n_d))
    f = 0._dp
    Allocate(f_list(1:2,1:n_E,1:n_d))
    f_list = 0._dp
    !Run n_p histories
    Do p = 1,n_p
        If (Mod(p,1000) .EQ. 0) Write(*,'(A,F6.2,A)',ADVANCE='NO') creturn,100._dp*Real(p,dp)/Real(n_p,dp),'%'
        ! start a photon at the source
        r = (/ 0._dp , 0._dp , 0._dp /)
        hv = hv0_pdf(RNG)
        omega_hat = Isotropic_Omega_Hat(RNG)
        t = 0._dp
        w = 1._dp
        r0 = 0._dp
        zeta0 = 1._dp
        Do
            ! move the photon to the location of the next interaction
            tau = -Log(1._dp - RNG%r())
            sig = interp_sig(hv,n_sig,hv_sigs,sigs)
            sigT = Sum(sig)
            d = tau / (rho0 * sigT)
            r = r + d * omega_hat
            t = t + d / speed_c
            !check for and tally crossings of the detector surfaces
            r1 = Vector_Length(r)
            Do s = 1,n_d
                If (r0 .GT. Rd(s)) Then !start outside...
                    Rca = r0 * Sqrt(1._dp - zeta0**2)
                    If (r1 .LT. Rd(s)) Then !finish inside, one crossing
                        zeta = Sqrt(1._dp - (Rca/Rd(s))**2)
                        If (p_listed(s)+1 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                        p_listed(s) = p_listed(s) + 1
                        p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                    Else If (zeta0 .LT. 0._dp) Then !finish outside, but initially traveling inward...
                        If (Rca .LT. Rd(s)) Then !radius of close approach is inside, potential for two crossings exists...
                            dca = Sqrt(r0**2 - Rca**2)
                            If (dca .LT. d) Then !distance traveled exceeds that required to reach Rca, two crossings
                                zeta = Sqrt(1._dp - (Rca/Rd(s))**2)
                                If (p_listed(s)+2 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                                p_listed(s) = p_listed(s) + 1
                                p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                                dz = Rd(s) - Rca
                                zeta = Sqrt(2._dp*Rca*dz + dz**2 ) / Rd(s)
                                p_list(2,p_listed(s),s) = p_list(2,p_listed(s),s) + w*Min(1./zeta,10._dp)
                            End If
                        End If
                    End If
                Else !start inside...
                    If (r1 .GT. Rd(s)) Then !finish outside, one crossing
                        If (zeta0 .GT. 0._dp) Then
                            dz = Rd(s) - r0
                            zeta = Sqrt((r0*zeta0)**2 + 2._dp*r0*dz + dz**2) / Rd(s)
                        Else
                            Rca = r0 * Sqrt(1._dp - zeta0**2)
                            dz = Rd(s) - Rca
                            zeta = Sqrt(2._dp*Rca*dz + dz**2) / Rd(s)
                        End If
                        If (p_listed(s)+1 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                        p_listed(s) = p_listed(s) + 1
                        p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                    End If
                End If
            End Do
            r0 = r1
            ! interact the photon in the medium
            xi = RNG%r() * sigT
            sig_sum = 0._dp
            Do interaction = 1,4
                sig_sum = sig_sum + sig(interaction)
                If (sig_sum .GT. xi) Exit
            End Do
            If (interaction .EQ. photoelectric_absorption) Exit
            Select Case (interaction)
                Case (incoherent_scatter,coherent_scatter)
                    If (interaction .EQ. incoherent_scatter) Then
                        alpha = hv / m0c2
                        mu0 = mu_KN(RNG,alpha)
                        hv = hv / (1._dp + alpha*(1._dp - mu0))
                    Else !coherent_scatter
                        mu0 = mu_Thomp(RNG)
                    End If
                    w0 = Isotropic_Azimuth(RNG)
                    If ( Abs(Dot_Product(omega_hat,Z_hat)) .LT. 0.5_dp ) Then
                        b_hat = Unit_Vector(Cross_Product(omega_hat,Z_hat))
                    Else
                        b_hat = Unit_Vector(Cross_Product(omega_hat,Y_hat))
                    End If
                    c_hat = Cross_Product(omega_hat,b_hat)
                    nu = Sqrt(1._dp - mu0**2)
                    omega_hat = mu0*omega_hat + nu*Cos(w0)*b_hat + nu*Sin(w0)*c_hat
                Case (pair_production)
                    omega_hat = Isotropic_Omega_hat(RNG)
                    hv = m0c2
                    w = w * 2._dp
            End Select
            zeta0 = Dot_Product(Unit_Vector(r),omega_hat)
        End Do
        !collect and record tallies from the history that just finished
        Do s = 1,n_d
            If (p_listed(s) .EQ. 0) Cycle
            i = p_listed(s)
            j = n_E
            fp = 0._dp
            Do
                j = Bisect_Search(p_list(1,i,s),E_list(0:j),n_E)
                Do
                    If (i .GT. 1) Then
                        If (p_list(1,i-1,s) .GT. E_list(j-1)) Then
                            p_list(1,i-1,s) = p_list(1,i-1,s) + p_list(1,i,s)
                            i = i - 1
                        Else
                            Exit
                        End If
                    Else
                        Exit
                    End If
                End Do
                fp = fp + p_list(2,i,s)
                f_list(:,j,s) = f_list(:,j,s) + (/ p_list(2,i,s) , p_list(2,i,s)**2 /)
                If (i .EQ. 1) Exit
                i = i - 1
            End Do
            f(:,s) = f(:,s) + (/ fp , fp**2 /)
        End Do
        p_listed = 0
        p_list = 0._dp
    End Do
    Write(*,*)
    Write(*,*)
    Do i = 1,n_d
        Write(*,'(2ES23.15)') f(1,i) / (n_p*4._dp*Pi*Rd(i)**2) , coeff_of_var(n_p,f(1,i),f(2,i))
    End Do
End Subroutine Transport_Photons

Subroutine Read_MAC(file_name,n,ln_hv,ln_sig)
    Use Kinds, Only: dp
    Use FileIO_Utilities, Only: Output_Message
    Implicit None
    Character(*), Intent(In) :: file_name
    Integer, Intent(Out) :: n
    Real(dp), Allocatable, Intent(Out) :: ln_hv(:)
    Real(dp), Allocatable, Intent(Out) :: ln_sig(:,:)
    Integer :: i
    Character(5) :: ctrash
    Real(dp) :: trash
    Integer :: unit, stat

    !open the file
    Open(NEWUNIT = unit , FILE = file_name , STATUS = 'OLD' , ACTION = 'READ' , IOSTAT = stat)
    If (stat .NE. 0) Call Output_Message('MACs: Read_MAC:  Error opening file '//file_name//', stat=',stat)
    !count the number of entries
    n = 0
    Do
        Read(unit,'(A5,2ES13.5E2)' , IOSTAT = stat) ctrash,trash,trash
        If (stat .LT. 0) Exit
        n = n + 1
    End Do
    Rewind(unit)
    !Allocate arrays
    Allocate(ln_hv(1:n))
    ln_hv = -1._dp
    Allocate(ln_sig(1:4,1:n))
    ln_sig = -1._dp
    !Read values
    Do i = 1,n
        Read(unit,'(7ES10.3E2)') ln_hv(i) ,    & !energy [MeV]
                               & ln_sig(3,i) , & !MAC for coherent scatter [g/cm^2]
                               & ln_sig(2,i) , & !MAC for incoherent scatter [g/cm^2]
                               & ln_sig(1,i) , & !MAC for photoelectric absorption [g/cm^2]
                               & ln_sig(4,i) , & !MAC for pair production (nuclear) [g/cm^2]
                               & trash           !MAC for pair production (atomic) [g/cm^2]
        ln_sig(4,i) = ln_sig(4,i) + trash  !merge pair production (nuclear/atomic) MACs
    End Do
    !close the file
    Close(unit)
    !convert to keV
    ln_hv = ln_hv * 1000._dp
    !convert to log scale
    ln_hv = Log(ln_hv)
    Where (ln_sig .NE. 0._dp)  !there are likely zeros in the pair production column
        ln_sig = Log(ln_sig)
    Elsewhere
        ln_sig = -33._dp
    End Where
End Subroutine

Function hv0_pdf_Line_N(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 7
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_N
Function hv0_pdf_Line_Fe(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 26
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_Fe
Function hv0_pdf_Line_Pb(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 82
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_Pb
Function hv0_pdf_Line_5(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG

    hv = 5.E3_dp
End Function hv0_pdf_Line_5

Function Planck(hv,kT) Result(p)
    !Returns the intensity at E [keV] of a normalized Plankian blackbody with temperature kT [keV]
    Use Kinds, Only: dp
    Use Global, Only: Pi
    Implicit None
    Real(dp) :: p
    Real(dp), Intent(In) :: hv
    Real(dp), Intent(In) :: kT
    Real(dp):: b,b_sq
    Real(dp), Parameter :: c1 =  1._dp / 47900160._dp
    Real(dp), Parameter :: c2 = -1._dp / 1209600._dp
    Real(dp), Parameter :: c3 =  1._dp / 30240._dp
    Real(dp), Parameter :: c4 = -1._dp / 720._dp
    Real(dp), Parameter :: c5 =  1._dp / 12._dp
    Real(dp), Parameter :: c6 = -1._dp / 2._dp

    b = hv / kT
    If (b .LE. 0._dp) Then
        p = 0._dp
        Return
    Else If (b .LT. 0.1_dp) Then
        !for small values of z use a series expansion approximation
        !No doubt, a more compact or recursive form exists, but the following correctly accumulates ~20 digits of precision
        b_sq = b**2
        p = b_sq * c1 + c2
        p = p * b_sq + c3
        p = p * b_sq + c4
        p = p * b_sq + c5
        p = p * b + c6
        p = p * b + 1._dp
        p = p * b_sq * 15._dp / (kT*(Pi**4))
    Else
        !Bridgman 5-7
        p = 15._dp * hv**3 / ((Exp(b) - 1._dp) * (kT*pi)**4)
    End If
End Function Planck

Function hv0_pdf_Planck_N(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 7
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_N
Function hv0_pdf_Planck_Fe(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 26
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_Fe
Function hv0_pdf_Planck_Pb(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 82
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_Pb

Function Isotropic_Azimuth(RNG) Result(w)
    !Returns w (or omega), an azimuthal scattering angle uniformly distributed [0,2pi)
    Use Kinds, Only: dp
    Use Global, Only: TwoPi
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: w
    Type(ACORN_Type), Intent(InOut) :: RNG

    w = TwoPi * RNG%r()
End Function Isotropic_Azimuth

Function Isotropic_mu(RNG) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta),
    ! uniformly distributed [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG

    mu = 2._dp * RNG%r() - 1._dp
End Function Isotropic_mu

Function mu_w_2_OmegaHat(xi,w) Result(Omegahat)
    !Converts angles (polar cosine and azimuthal rotation) to a unit vector
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: OmegaHat(1:3)
    Real(dp), Intent(In) :: xi !cosine of the polar scattering angle
    Real(dp), Intent(In) :: w  !azimuthal scattering angle
    Real(dp) :: nu,mu,eta

    nu = Sqrt(1._dp - xi**2)
    mu = nu * Cos(w)
    eta = nu * Sin(w)
    OmegaHat = (/ mu, eta, xi /)
End Function mu_w_2_OmegaHat

Function Isotropic_Omega_hat(RNG) Result(OmegaHat)
    !Reutrns an isotropic random direction
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp):: OmegaHat(1:3)
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp):: xi  ! z direction cosine = cos(theta), theta is angle from zHat
    Real(dp):: w   ! angle from xHat toward yHat around zHat

    xi = Isotropic_mu(RNG)
    w = Isotropic_Azimuth(RNG)
    OmegaHat = mu_w_2_OmegaHat(xi,w)
End Function Isotropic_Omega_hat

Function Bisect_Search(x, xx, n) Result (j)
    ! Finds index j, such that x is in the interval ( xx(j-1),xx(j) ]
    ! The list xx must be sorted in acsending order
    ! j=1 or j=n+1 indicates that x is not in the range [ xx(1),xx(n) ]
    Use Kinds, Only: dp
    Implicit None
    Integer :: j
    Real(dp), Intent(In) :: x
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: xx(1:n)
    Integer :: jl,jh,jm

    !First, check for endpoints
    If (x .EQ. xx(1)) Then
        j = 2
        Return
    Else If (x .EQ. xx(n)) Then
        j = n
        Return
    End If
    !Otherwise, search the list
    jl = 0
    jh = n + 1
    Do
        If (jh - jl .LE. 1) Exit
        jm = (jh + jl) / 2  !integer arithmetic
        If (x .LT. xx(jm)) Then
            jh = jm
        Else
            jl = jm
        End If
    End Do
    j = jh
End Function Bisect_Search

Elemental Function Lin_Interp(x,x1,x2,y1,y2) Result(y)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: y
    Real(dp), Intent(In) :: x
    Real(dp), Intent(In) :: x1,x2
    Real(dp), Intent(In) :: y1,y2

    y = y1 + (y2-y1) * (x-x1) / (x2-x1)
End Function Lin_Interp

Function interp_sig(hv,n_sig,hv_sig,sig)
    Use Kinds, Only: dp
    Use Global, Only: m0c2
    Implicit None
    Real(dp) :: interp_sig(1:4)
    Real(dp), Intent(In) :: hv
    Integer, Intent(In) :: n_sig
    Real(dp), Intent(In) :: hv_sig(1:n_sig)
    Real(dp), Intent(In) :: sig(1:4,1:n_sig)
    Real(dp), Parameter :: pp_thresh = 2._dp * m0c2
    Real(dp) :: ln_hv
    Integer :: j

    ln_hv = Log(hv)
    j = Bisect_Search(ln_hv,hv_sig,n_sig)
    If (j .EQ. 0) j = 1
    interp_sig = Exp( Lin_Interp(ln_hv,hv_sig(j-1),hv_sig(j),sig(:,j-1),sig(:,j)) )
    If (hv .LT. pp_thresh) interp_sig(4) = 0._dp
End Function interp_sig

Function mu_Thomp(RNG) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta),
    ! distributed on the Thompson distribution on [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp) :: xi
    Real(dp) :: a
    Real(dp), Parameter :: one_third = 1._dp / 3._dp
    Real(dp), Parameter :: two_thirds = 2._dp / 3._dp

    xi = RNG%r()
    a = 2._dp - 4._dp * xi + Sqrt(5._dp + 16._dp*xi*(xi-1._dp))
    mu = (1._dp - a**two_thirds) / (a**one_third)
End Function mu_Thomp

Function mu_KN(RNG,alpha) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta),
    ! distributed on the Klein-Nishina Distribution on [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp), Intent(In) :: alpha ![] energy of the incident photon, in rest mass energy of an electron
    Real(dp) :: kn  !Klein Nishina distribution evaluated at a proposed mu
    Real(dp) :: mu2p1 !(mu**2 + 1._dp) intermediate value for computing kn
    Real(dp) :: amamup1 !( 1._dp + a * (1._dp - mu) ) intermediate value for computing kn

    Do  !sampling by geometric rejection
        !sample mu from Thompson scattering distribution
        mu = mu_Thomp(RNG)
        !compute the value of the KN distribution at this mu
        !This is the KN distributon scaled on 0 < f(mu) < TWO  (using this scaled value eliminates a multiplication here and a
        !                                                       multiplication in the following check for acceptance)
        mu2p1 = 1._dp + mu**2
        amamup1 = 1._dp + alpha * (1._dp - mu)
        kn = mu2p1 * ( 1._dp + (alpha * (1._dp - mu))**2 / (mu2p1 * amamup1) ) / amamup1**2
        !If the probability of KN at this mu exceeds the probability of the sampled dist, then accept the value and exit
        If ( kn .GT. mu2p1*RNG%r() ) Exit
        !Otherwise, cycle and repeat the sampling process
    End Do
End Function mu_KN

Pure Function Vector_Length(v) Result(d)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: d
    Real(dp), Intent(In) :: v(1:3)

    d = Sqrt(Sum(v*v))
End Function Vector_Length

Pure Function Unit_Vector(v) Result(hat)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: hat(1:3)
    Real(dp), Intent(In) :: v(1:3)

    hat = v / Vector_Length(v)
End Function Unit_Vector

Pure Function Cross_Product(v,w) Result(u)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: u(1:3)
    Real(dp), Intent(In) :: v(1:3),w(1:3)

    u = (/ v(2)*w(3) - v(3)*w(2) , &
         & v(3)*w(1) - v(1)*w(3) , &
         & v(1)*w(2) - v(2)*w(1) /)
End Function Cross_Product

Subroutine Resize_p_list(nd,n,p_list)
    Use Kinds, Only: dp
    Implicit None
    Integer, Intent(In) :: nd
    Integer, Intent(InOut) :: n
    Real(dp), Allocatable, Intent(InOut) :: p_list(:,:,:)
    Real(dp), Allocatable :: p_swap(:,:,:)

    Allocate(p_swap(1:2,1:n,1:nd))
    p_swap = p_list
    Deallocate(p_list)
    Allocate(p_list(1:2,1:2*n,1:nd))
    p_list = 0._dp
    p_list(1:2,1:n,1:nd) = p_swap
    Deallocate(p_swap)
    n = 2 * n
End Subroutine Resize_p_list

Function std_err( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: std_err
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    If (n .EQ. 1) Then
        std_err = mu
    Else
        std_err = Sqrt( (mu_sq - mu**2 / Real(n,dp)) / Real(n-1,dp) )
    End If
End Function std_err

Function std_err_otm( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: std_err_otm
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    std_err_otm = std_err(n,mu,mu_sq) / Sqrt(Real(n,dp))
End Function std_err_otm

Function coeff_of_var( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: coeff_of_var
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    If (n .EQ. 1) Then
        coeff_of_var = 1._dp
    Else
        coeff_of_var = std_err_otm(n,mu,mu_sq) / (mu/Real(n,dp))
    End If
End Function coeff_of_var

End Program

max_stars_repo_name: wtdailey/XrayToolKit

max_issues_repo_name: wtdailey/XrayToolKit

max_forks_repo_name: wtdailey/XrayToolKit


text: #redirect Hoa Viet

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


SUBROUTINE classify_by_eigenvalues(classification,classnum, values,threshold)
!--------------------------------------------------------------------------
! This subroutine compares three eigenvalues to a threshold, and determines
! how many eigenvalues exceed that threshold (N_above +1)
!
! Classification = N_above +1
! (+1 because fortran can't deal with zero array indices)
!---------------------------------------------------------------------------
implicit none

integer,intent(inout) :: classification
real, intent(in) :: threshold
integer,dimension(4), intent(inout) :: classnum
real,dimension(3) :: values

integer :: k

classification = 1

! Find classification
do k=1,3
   if(values(k)<threshold) then
      classification = classification +1
   endif
enddo

! Add to totals for each class
classnum(classification) = classnum(classification)+1

return

END SUBROUTINE classify_by_eigenvalues


!-----------------------------------------------------------
! Subroutine to write radial dumps
!-----------------------------------------------------------

subroutine write_dump

  use gravdata
  use magdata
  use planetdata
  use unitdata
  use winddata

  implicit none

  real(kind=8) :: sig_max,mdisk, tot_lumin,tcoolmin
  real :: mdot_grav, mdot_mag, grav_max,mag_max, mmag,mgrav
  integer :: i,ifirst, iplanet

  write (*,103) 'Output at time t = ',t/yr, '  years, dt = ',dt/yr, ' years'
102 format (11(1PE15.5))
103 format (A,1PE11.3,A,1PE11.3,A)
104 format (4(1PE15.5))
110 format(10(1PE15.5))
111 format (11(1PE15.5))

  ! Calculate disc properties and spectrum

  call disc_properties
  call luminosity(tot_lumin,spectrum,Tc,tau)

  write(fileno, snapshotformat)snapshotcounter

  fileno = TRIM(fileno)
! write out disc profiles
  OPEN(iprof, file=TRIM(prefix)//'_profile.'//fileno, status='unknown')
  write(iprof,*) t/yr, nrgrid
  do i = isr, ier
     write (iprof,102) rz(i)/AU,sigma(i), cs(i), kappa(i), gamma(i),mu(i),Tc(i), &
    tau(i), nu_tc(i), alpha_g(i),Q(i)
  enddo

  close(iprof)

  OPEN(ispec, file=TRIM(prefix)//'_spectrum.'//fileno, status='unknown')
  write(ispec,*) 1000
  do i = 1, 1000
     write(ispec,*) spectrum(1,i), spectrum(2,i)
  enddo

  close(ispec)

if (layerchoice=='y') then
OPEN(iprof,file=TRIM(prefix)//'layer.'//fileno,status='unknown')
WRITE(iprof,*) t/yr,nrgrid
do i=isr,ier
write(iprof,111) rz(i)/AU, sigma(i), sigma_m(i), sigma_tot(i), cs_m(i),&
    kappa_m(i),gamma_m(i), mu_m(i),tau_m(i),nu_m(i),alpha_m
enddo
close(iprof)
endif

if(planetchoice=='y') then
open(iprof, file=TRIM(prefix)//'_planets.'//fileno,status='unknown')
write(iprof,*)t/yr,nplanet,nactive
do iplanet=1,nplanet

    write(iprof,*) alive(iplanet),mp(iplanet)/mjup, &
         ap(iplanet)/AU, tmig(iplanet)/yr
enddo
close(iprof)

open(75,file=TRIM(prefix)//'_torque.'//fileno,status='unknown')
do i=isr,ier
   write(75,*) rz(i)/AU, 0.1*H(i)/rz(i), lambdaI(1,i), lambdaII(1,i), &
        total_planet_torque(i), torque_term(i), total_planet_torque(i+1)-total_planet_torque(i)
enddo
close(75)

endif

! Compute disk mass and maximum surface density
! Also compute radially averaged accretion rates

mdisk = 0.0d0
mgrav= 0.0d0
mmag = 0.0d0

grav_max = 0.0d0
mag_max = 0.0d0
sig_max = 0.0d0

do i = isr, ier
mmag = mmag + 2.0d0*pi*rz(i)*(rf(i+1)-rf(i))*sigma_m(i)
mgrav = mgrav + 2.0d0*pi*rz(i)*(rf(i+1)-rf(i))*sigma(i)
mdisk = mdisk+2.0d0*pi*rz(i)*(rf(i+1)-rf(i))*sigma_tot(i)

grav_max = max(grav_max,sigma(i))
mag_max = max(mag_max,sigma_m(i))
sig_max = max(sig_max,sigma_tot(i))

mdot_grav = mdot_grav + 3.0*pi*nu_tc(i)*sigma(i)/REAL(nrgrid)
mdot_mag = mdot_mag + 3.0*pi*nu_m(i)*sigma_m(i)/REAL(nrgrid)

mdot_wind = mdot_wind + sigdot_wind(i)*twopi*rz(i)/drzm1(i)

enddo

mdisk = mdisk/solarmass
mgrav = mgrav/solarmass
mmag = mmag/solarmass

mdot_grav = mdot_grav/msolyr
mdot_mag = mdot_mag/msolyr
mdot_wind = mdot_wind/msolyr

! Write out snapshot data (disk mass, sigma, total luminosity)
write(itime,111) t/yr, dt/yr, mdisk, tot_lumin, sig_max,mgrav,mmag,grav_max,mag_max,mdot_grav,mdot_mag,mdot_wind
call flush(itime)

  return

end  subroutine write_dump


text:
c     Calculate determinant of square matrix
c     Author: Pearu Peterson, March 2002
c
c     prefixes: d,z,s,c   (double,complex double,float,complex float)
c     suffixes: _c,_r     (column major order,row major order)

      subroutine ddet_c(det,a,n,piv,info)
      integer n,piv(n),i
      double precision det,a(n,n)
cf2py intent(in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument double*,double*,int*,int*,int*
      external dgetrf
      call dgetrf(n,n,a,n,piv,info)
      det = 0d0
      if (info.ne.0) then
         return
      endif
      det = 1d0
      do 10,i=1,n
         if (piv(i).ne.i) then
            det = -det * a(i,i)
         else
            det = det * a(i,i)
         endif
 10   continue
      end

      subroutine ddet_r(det,a,n,piv,info)
      integer n,piv(n)
      double precision det,a(n,n)
cf2py intent(c,in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument double*,double*,int*,int*,int*
      external ddet_c
      call ddet_c(det,a,n,piv,info)
      end

      subroutine sdet_c(det,a,n,piv,info)
      integer n,piv(n),i
      real det,a(n,n)
cf2py intent(in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument float*,float*,int*,int*,int*
      external sgetrf
      call sgetrf(n,n,a,n,piv,info)
      det = 0e0
      if (info.ne.0) then
         return
      endif
      det = 1e0
      do 10,i=1,n
         if (piv(i).ne.i) then
            det = -det * a(i,i)
         else
            det = det * a(i,i)
         endif
 10   continue
      end

      subroutine sdet_r(det,a,n,piv,info)
      integer n,piv(n)
      real det,a(n,n)
cf2py intent(c,in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument float*,float*,int*,int*,int*
      external sdet_c
      call sdet_c(det,a,n,piv,info)
      end

      subroutine zdet_c(det,a,n,piv,info)
      integer n,piv(n),i
      complex*16 det,a(n,n)
cf2py intent(in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument complex_double*,complex_double*,int*,int*,int*
      external zgetrf
      call zgetrf(n,n,a,n,piv,info)
      det = (0d0,0d0)
      if (info.ne.0) then
         return
      endif
      det = (1d0,0d0)
      do 10,i=1,n
         if (piv(i).ne.i) then
            det = -det * a(i,i)
         else
            det = det * a(i,i)
         endif
 10   continue
      end

      subroutine zdet_r(det,a,n,piv,info)
      integer n,piv(n)
      complex*16 det,a(n,n)
cf2py intent(c,in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument complex_double*,complex_double*,int*,int*,int*
      external zdet_c
      call zdet_c(det,a,n,piv,info)
      end

      subroutine cdet_c(det,a,n,piv,info)
      integer n,piv(n),i
      complex det,a(n,n)
cf2py intent(in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument complex_float*,complex_float*,int*,int*,int*
      external cgetrf
      call cgetrf(n,n,a,n,piv,info)
      det = (0e0,0e0)
      if (info.ne.0) then
         return
      endif
      det = (1e0,0e0)
      do 10,i=1,n
         if (piv(i).ne.i) then
            det = -det * a(i,i)
         else
            det = det * a(i,i)
         endif
 10   continue
      end

      subroutine cdet_r(det,a,n,piv,info)
      integer n,piv(n)
      complex det,a(n,n)
cf2py intent(c,in,copy) :: a
cf2py intent(out) :: det,info
cf2py integer intent(hide,cache),depend(n),dimension(n) :: piv
cf2py integer intent(hide),depend(a) :: n = shape(a,0)
cf2py check(shape(a,0)==shape(a,1)) :: a
cf2py callprotoargument complex_float*,complex_float*,int*,int*,int*
      external cdet_c
      call cdet_c(det,a,n,piv,info)
      end

max_stars_repo_name: lesserwhirls/scipy-cwt

max_issues_repo_name: lesserwhirls/scipy-cwt

max_forks_repo_name: lesserwhirls/scipy-cwt


! for more info refer : https://pages.mtu.edu/~shene/COURSES/cs201/NOTES/chap02/funct.html
program intrinsicfunction

! CTM is most of the Mathematical function the argument must be real

print *, "abs(-5)            : ", abs(-5); ! argument can be integer/real
print *, "abs(25.14)         : ", abs(25.14);
print *, "mod(10, 2)         : ", mod(10, 2);
print *, "mod(17, 3)         : ", mod(10, 2);
print *, "sqrt(16.0)         : ", sqrt(16.0);
print *, "sqrt(13)           : ", sqrt(13.5);
print *, "exp(3)             : ", exp(3.0);    ! base e (euler e = 2.718)
print *, "log(10)            : ", log(10.0);
print *, "int(12.41)         : ", int(12.41); ! gives integer part
print *, "nint(11.43)        : ", nint(11.43);
print *, "nint(11.63)        : ", nint(11.63);
print *, "fraction(11.43)    : ", fraction(11.43);
print *, "floor(11.43)       : ", floor(11.43);
print *, "floor(11.93)       : ", floor(11.93);
print *, "real(10)           : ", real(10);
print *, "max(3, 44, 54)     : ", max(3, 44, 54);
print *, "min(3, 44, 54)     : ", min(3, 44, 54);
print *, "max(3.4, 6.1, 5.6) : ", max(3.4, 6.1, 5.6);
print *, "min(3.4, 6.1, 5.6) : ", min(3.4, 6.1, 5.6);

end program intrinsicfunction


MODULE Meltingmod
    IMPLICIT NONE
CONTAINS

    SUBROUTINE model_melting(jul, &
        Tmf, &
        DKmax, &
        Kmin, &
        tavg, &
        M)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: jul
        REAL, INTENT(IN) :: Tmf
        REAL, INTENT(IN) :: DKmax
        REAL, INTENT(IN) :: Kmin
        REAL, INTENT(IN) :: tavg
        REAL, INTENT(OUT) :: M
        REAL:: K
        !- Name: Melting -Version: 1.0, -Time step: 1
        !- Description:
    !            * Title: snow in the process of melting
    !            * Author: STICS
    !            * Reference: doi:http://dx.doi.org/10.1016/j.agrformet.2014.05.002
    !            * Institution: INRA
    !            * Abstract: It simulates snow in the process of melting
        !- inputs:
    !            * name: jul
    !                          ** description : current day of year for the calculation
    !                          ** inputtype : variable
    !                          ** variablecategory : auxiliary
    !                          ** datatype : INT
    !                          ** default : 0
    !                          ** min : 0
    !                          ** max : 366
    !                          ** unit : d
    !                          ** uri :
    !            * name: Tmf
    !                          ** description : threshold temperature for snow melting
    !                          ** inputtype : parameter
    !                          ** parametercategory : constant
    !                          ** datatype : DOUBLE
    !                          ** default : 0.5
    !                          ** min : 0.0
    !                          ** max : 1.0
    !                          ** unit : degC
    !                          ** uri :
    !            * name: DKmax
    !                          ** description : difference between the maximum and the minimum melting rates
    !                          ** inputtype : parameter
    !                          ** parametercategory : constant
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 5000.0
    !                          ** unit : mmW/degC/d
    !                          ** uri :
    !            * name: Kmin
    !                          ** description : minimum melting rate on 21 December
    !                          ** inputtype : parameter
    !                          ** parametercategory : constant
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 5000.0
    !                          ** unit : mmW/degC/d
    !                          ** uri :
    !            * name: tavg
    !                          ** description : average temperature
    !                          ** inputtype : variable
    !                          ** variablecategory : auxiliary
    !                          ** datatype : DOUBLE
    !                          ** default : 0.0
    !                          ** min : 0.0
    !                          ** max : 100.0
    !                          ** unit : degC
    !                          ** uri :
        !- outputs:
    !            * name: M
    !                          ** description : snow in the process of melting
    !                          ** variablecategory : rate
    !                          ** datatype : DOUBLE
    !                          ** min : 0.0
    !                          ** max : 500.0
    !                          ** unit : mmW/d
    !                          ** uri :
        M = 0.0
        K = DKmax / 2.0 * (-SIN((2.0 * 3.14159265 * REAL(jul) / 366.0 + (9.0  &
                / 16.0 * 3.14159265)))) + Kmin + (DKmax / 2.0)
        IF(tavg .GT. Tmf) THEN
            M = K * (tavg - Tmf)
        END IF
    END SUBROUTINE model_melting

END MODULE