!  Paralleltest1.f90
!
! Here are my results, comparing Fortran without OpenMP (i e, single thread), Fortran with OpenMP, and C++ using my own thread control (no OpenMP or other library).
!
! Fortran,  Release:    8 threads 10000 radii  Elapsed=   20.0340 sec  CPU time=   144.3945 sec sum= 2.6670250794563941E+02
!                       1 threads 10000 radii  Elapsed=    5.7480 sec  CPU time=     5.6940 sec sum= 2.6670250794563941E+02
!
! C++,      Release:    8 threads 10000 radii  Elapsed=    5.7300 sec  CPU time=        n/a sec sum= 2.6670250794564652e+02
!
! Compiled under Visual Studio as x64 release build, with option /Qopenmp.
!
! 3.4 GHz Haswell (8 logical CPUs), Windows 7 64-bit, Visual Studio 2017 15.6.7, Intel XE2018 update 1.
!
!
  MODULE ParallelTestData
    INTEGER :: nRadii, nSpecies
    real*8, DIMENSION(3,10000) :: A1, A2, B
    real*8, DIMENSION(10000,3) :: U, V, W
    real*8, DIMENSION(0:10000+1) :: Cu, Cv, Cw, Ca
    real*8 :: dt, omSqRun
    real*8, DIMENSION(0:10001,28) :: Cm
    LOGICAL :: nonIdeal
    real*8, DIMENSION(28) :: diff0, sed0, BM1, Ks
    INTEGER, DIMENSION(0:3, 32) :: kChg
  END MODULE ParallelTestData

  PROGRAM Paralleltest1
    USE ParallelTestData
    USE omp_lib

    IMPLICIT NONE
    
    real*8, EXTERNAL :: DiffCoeff, SedCoeff, ZCalc

    INTEGER :: iRun, iRadius, iSpecies, nThreads, nRSpacings, iDt, iThread, ii
    real*8, DIMENSION(0:10000+1) :: Dradius,sRadius
    real*8, DIMENSION(10000)     :: Z
    real*8, DIMENSION(3,10000+1) :: concTotal
    real*8, DIMENSION(3,10000)   :: G
    real*8, DIMENSION(3)         :: rSpacings

    INTEGER :: iDt_loop, iDt_num, iDt_one
    INTEGER :: countStart,count_rate,countEnd, op_n
    real*8  :: startTime, endTime, CPUTime, elsapsedTime, sum_dot, sumx
    real*8  :: times(2,5)
    CHARACTER(LEN=100) :: msg
    logical :: no_hyper_threading = .true.

    integer*8 :: count
!
    nThreads = -1
!
    iDt_loop = 100000
    nRadii   = 10000
    nSPecies = 1
    iSpecies = 1
    iRun     = 1
    diff0(1) = 1e-13
    sed0(1)  = 1e-7
    nRSpacings=1; rSpacings(1)=1.705668493358788D-003
    kChg(0,1)=1; kChg(1,1)=nRadii
    Cm(1:nRadii, 1)=1.23
    nonIdeal=.FALSE.

    omSqRun = 1

    CALL Grid(nRSpacings, rSpacings, iSpecies, iRun)
    CALL CUVWCalc(iSpecies, iRun, concTotal)

    Dradius(0)=0; Dradius(nRadii+1)=0
    sRadius(0)=0; sRadius(nRadii+1)=0
    dt=0.01/1.02

      times = 0
      call delta_step (times(1,1))    ! initialise timer
      times = 0
      
 CALL SYSTEM_CLOCK (countStart,count_rate)
 CALL CPU_TIME     (startTime)

    sum_dot = 0
    sumx    = 0
    count   = 0
    iDt_num = 0
    iDt_one = 0

    if ( no_hyper_threading ) then
       nThreads = omp_get_max_threads ()
       if ( nThreads <  8 ) nThreads = nThreads-1    ! spare thread for 4-thread
       if ( nThreads >= 8 ) nThreads = nThreads/2    ! assume hyperthreading so n/2
       write (*,*) ' Threads reduced to',nThreads
       call omp_set_num_threads (nThreads)
    end if

!$OMP PARALLEL DEFAULT (SHARED)               &
!$OMP  PRIVATE (iThread,iRadius,ii)           &
!$OMP  REDUCTION (+ : sumx)                   &
!$OMP  REDUCTION (+ : count)                  &
!$OMP  REDUCTION (+ : iDt_num)                &
!$OMP  REDUCTION (+ : op_n)                   

!    OMP  LASTPRIVATE (nThreads)                 &

    iThread  = omp_get_thread_num  ()

!    OMP  ATOMIC
!$OMP  SINGLE
    nThreads = omp_get_num_threads ()
!$OMP  END SINGLE        
    
    DO iDt=1,iDt_loop
        iDt_num = iDt_num + 1

!$OMP  SINGLE
        dt=dt*1.00002
        iDt_one = iDt_one + 1
!$OMP  END SINGLE        
!
!$OMP  DO schedule(static)
        L_R1_Parallel: &
          DO iRadius = 1, nRadii
            Dradius(iRadius) = DiffCoeff (iSpecies, iRadius, concTotal(1:3, iRadius)) ! diffusion coefficient at iRadius
            sRadius(iRadius) = SedCoeff  (iSpecies, iRadius, concTotal(1:3, iRadius)) ! if compression, this will account for it
            ii = iRadius
          END DO L_R1_Parallel
!$OMP  END DO

      if ( iDt==iDt_loop ) write (*,*) iThread,NThreads,ii,iRadius

      IF ( iThread==0) call delta_step (times(1,1))

!$OMP  DO schedule(static)
        L_R2_Parallel: &
            DO iRadius = 1, nRadii
                Z(iRadius)   = ZCalc (iRadius, iSpecies)
                G(1, iRadius)= Dradius(iRadius-1)*dt*A1(1, iRadius) + B(1,iRadius) - sRadius(iRadius-1)*omSqRun*dt*A2(1, iRadius)
                G(2, iRadius)= Dradius(iRadius)  *dt*A1(2, iRadius) + B(2,iRadius) - sRadius(iRadius)  *omSqRun*dt*A2(2, iRadius)
                G(3, iRadius)= Dradius(iRadius+1)*dt*A1(3, iRadius) + B(3,iRadius) - sRadius(iRadius+1)*omSqRun*dt*A2(3, iRadius)
                sumx         = sumx + Z(iRadius) * G(1, iRadius) + G(2, iRadius)* G(3, iRadius)
            END DO L_R2_Parallel
!$OMP  END DO
!
      IF ( iThread==0) call delta_step (times(1,2))

!$OMP SINGLE
      sum_dot = sum_dot                                       &
          + DOT_PRODUCT(Z(1:nRadii),   G(1,1:nRadii))         &
          + DOT_PRODUCT(G(2,1:nRadii), G(3,1:nRadii))
!$OMP END SINGLE

      IF ( iThread==0) call delta_step (times(1,3))

      op_n = 0
!$OMP DO  schedule (static)            
      DO iRadius = 1, nRadii
        op_n = op_n+1
      END DO
!$OMP END DO
!
      IF ( iThread==0) call delta_step (times(1,4))
      count = count+1

      IF ( iThread==0) call delta_step (times(1,5))
    END DO
!$OMP END PARALLEL
    
    CALL SYSTEM_CLOCK(countEnd)
    CALL CPU_TIME    (endTime)

    write (*,*) 'nThread =', nThreads
    write (*,*) 'iDt_num =', iDt_num
    write (*,*) 'iDt_one =', iDt_one
    write (*,*) 'count   =', count
    write (*,*) 'op_n    =', op_n
    write (*,*) 'sums    =', sum_dot
    write (*,*) 'sumx    =', sumx
    
    elsapsedTime = DBLE(countEnd-countStart)/count_rate
    CPUTime      = endTime-startTime
    WRITE(msg,"(I0,' threads ',I0,' radii  Elapsed=',F10.4,' sec  CPU time= ',F10.4,' sec sum=', ES23.16)") &
       nThreads, nRadii, elsapsedTime, CPUTime, sum_dot
    WRITE(*,"(A)") msg
!    
    nThreads=omp_get_num_threads()
    write (*,*) 'sums =', sum_dot, sumx, nThreads
    !READ(*,*)
    
    write (*,*) ' '
    write (*,*) '!$omp entries',iDt_loop
    write (*,*) 'omp loop size',op_n
    write (*,11) 'loop1', times(:,1), times(1,1)/dble(iDt_loop), '  function calls'
    write (*,11) 'loop2', times(:,2), times(1,2)/dble(iDt_loop), '  calculation'
    write (*,11) 'sum  ', times(:,3), times(1,3)/dble(iDt_loop), '  dot_product'
    write (*,11) 'loop3', times(:,4), times(1,4)/dble(iDt_loop), '  OMP overhead'
    write (*,11) 'timer', times(:,5), times(1,5)/dble(iDt_loop), '  timer overhead'
    write (*,11) 'Total', sum(times(1,1:5)), sum(times(2,1:5)), sum(times(1,1:5))/dble(iDt_loop), '  Total'
 11 format (a,3en12.2,a)
 
    END PROGRAM Paralleltest1

    FUNCTION ZCalc(iRadius, iSpecies)
    USE ParallelTestData
    IMPLICIT NONE
    real*8 :: ZCalc
    INTEGER, INTENT(IN) :: iRadius, iSpecies
    real*8 :: uCu, vCv, wCw, aCa, concDependent, nonIdealCorr


    ZCalc = B(1,iRadius)*Cm(iRadius-1, iSpecies)     &
          + B(2,iRadius)*Cm(iRadius,   iSpecies)     &
          + B(3,iRadius)*Cm(iRadius+1, iSpecies)

    IF(nonIdeal) THEN
        ! Now, the concentration-dependent part for diffusion
        !v6.51: called once        CALL CUVWCalc(iSpecies, concTotal, Cu, Cv, Cw, Ca)
        uCu = U(iRadius, 1)*Cu(iRadius-1) &
            + U(iRadius, 2)*Cu(iRadius)   &
            + U(iRadius, 3)*Cu(iRadius+1)
        vCv = V(iRadius, 1)*Cv(iRadius-1) &
            + V(iRadius, 2)*Cv(iRadius)   &
            + V(iRadius, 3)*Cv(iRadius+1)
        wCw = W(iRadius, 1)*Cw(iRadius-1) &
            + W(iRadius, 2)*Cw(iRadius)   &
            + W(iRadius, 3)*Cw(iRadius+1)
        concDependent = uCu + vCv + wCw ! this part is for D only
        aCa = A2(1, iRadius)*Ca(iRadius-1) & ! this part is for s
            + A2(2, iRadius)*Ca(iRadius)   &
            + A2(3, iRadius)*Ca(iRadius+1)

        !nonIdealCorr = dt*gammaC(iSpecies)*concDependent ! gammaC = diff0(iSpecies)*(-2*BM1(iSpecies))
        nonIdealCorr = dt*(diff0(iSpecies)*concDependent - sed0(iSpecies)*omSqRun*aCa)
        ZCalc = ZCalc + nonIdealCorr
    END IF ! nonIdeal
    RETURN
    END FUNCTION ZCalc

    ! DiffCoeff.f90: Calculates the concentration-independent or concentration-dependent diffusion coefficient

    FUNCTION DiffCoeff(iSpecies, iRadius, cTot)
    USE ParallelTestData
    IMPLICIT NONE

    real*8 :: DiffCoeff
    INTEGER, INTENT(IN) :: iSpecies, iRadius
    real*8, DIMENSION(3), INTENT(IN) :: cTot
    real*8 :: totNonIdeal

    DiffCoeff=diff0(iSpecies) ! constant independent of radius for a single species
    totNonIdeal=cTot(1)

    RETURN
    END

    FUNCTION SedCoeff(iSpecies, iRadius, cTot)

    USE ParallelTestData
    IMPLICIT NONE

    real*8 :: SedCoeff
    INTEGER, INTENT(IN) :: iSpecies, iRadius
    real*8, DIMENSION(3), INTENT(INOUT) :: cTot

    real*8 :: totNonIdeal
    INTEGER :: n1=0

    SedCoeff=sed0(iSpecies) ! constant independent of radius for a single species
    totNonIdeal=cTot(1)

    RETURN
    END

    SUBROUTINE CUVWCalc(iSpecies, iRun, concTotal)
    USE ParallelTestData
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: iSpecies, iRun
    real*8, DIMENSION(3, 10000+1), INTENT(IN) :: concTotal
    !real*8 :: sumSWeighted, sumDWeighted, sumC, sumNC, s0, D0, cN, niThermo, niHydro, totNonIdeal, sNonIdeal, DNonIdeal

    INTEGER :: iRadius, iRad, jSpecies, iList, iS !, messageBoxResult
    real*8, DIMENSION(1:10000) :: P, Q
    real*8, DIMENSION(0:10000+1) :: eta, eps ! the functions of concentration in s=s_0(1-ks*eps) and D=D0(1-kD*eta)

    DO iRadius=1,nRadii
        P(iRadius) = 2.*DOT_PRODUCT(BM1(1:nSpecies),Cm(iRadius,1:nSpecies))
    END DO


    DO iRadius=1,nRadii
        Q(iRadius) = DOT_PRODUCT(Ks(1:nSpecies),Cm(iRadius,1:nSpecies))
    END DO

    ! eta = 1 - D/D0 (T&H 1981 notation)
    eta(0)=0
    eta(1:nRadii) = (Q(1:nRadii) - P(1:nRadii)) / (1 + Q(1:nRadii))


    eta(nRadii+1)=0 ! extra 0-elements make the end cases 0 where they don't exist

    ! eps = 1 - s/s0  (T&H 1981 notation)
    eps(0)=0
    eps(1:nRadii) = Q(1:nRadii) / (1 + Q(1:nRadii))
    eps(nRadii+1)=0

    ! Debugging output if wanted

    DO iRadius=1,nRadii
        Cu(iRadius)=Cm(iRadius+1, iSpecies)*eta(iRadius) ! Todd & Haschemeyer
        Cv(iRadius)=Cm(iRadius,   iSpecies)*eta(iRadius)
        Cw(iRadius)=Cm(iRadius,   iSpecies)*eta(iRadius+1)

        Ca(iRadius)=Cm(iRadius, iSpecies)*eps(iRadius)
    END DO
    Cu(0)=0; Cu(nRadii+1)=0
    Cv(0)=0; Cv(nRadii+1)=0
    Cw(0)=0; Cw(nRadii+1)=0
    Ca(0)=0; Ca(nRadii+1)=0

    RETURN
    END SUBROUTINE CUVWCalc

    SUBROUTINE Grid(nRSpacings, rSpacings, iSpecies, iRun)

    USE ParallelTestData
    USE, INTRINSIC :: IEEE_Arithmetic

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: nRSpacings, iSpecies, iRun
    real*8, DIMENSION(3), INTENT(IN) :: rSpacings

    real*8, DIMENSION(10000) :: gradient, dRa
    real*8 :: dx0, dxMin, dx1, dx2, sumRGradients, a, r0, r1, rx, gMax, RCheck, gradientMin
    INTEGER :: iSpacing, iRadius, iRadius0, nGradZero, n, nGradLarge, debugOutputUnit, nGradZero1, nGradLarge1, nTotal, nNormal

    iRadius0=1; dx2=rSpacings(1); iSpacing=1

    dx1=dx2
    r0=5.9; r1=7.2
    nGradZero1=0; nGradLarge1=0

    L_R:DO iRadius=1, nRadii-1
        IF(iRadius>=kChg(iSpacing, iRun)) THEN
            iSpacing=iSpacing+1
            dx1=dx2
            dx2=rSpacings(iSpacing)
        END IF

        rx=r0+(iRadius-1)*dx2
        B(2, iRadius)  =(dx1+dx2)*(rx/3.-(dx1-dx2)/12.)
        B(3, iRadius)  =dx2*(rx+dx2/2.0)/6.0
        B(1, iRadius+1)=B(3, iRadius)

        A1(2, iRadius)  =(rx-dx1/2.)/dx1 + (rx+dx2/2.)/dx2
        A1(3, iRadius)  =-1.0*(rx+dx2/2.0)/dx2
        A1(1, iRadius+1)=A1(3, iRadius)

        A2(2, iRadius)  =((dx1*dx1-dx2*dx2)/4. - rx*(dx1+dx2))/3.0
        A2(3, iRadius)  =-1.0*(rx*rx/2.0+2.0*rx*dx2/3.0+dx2*dx2/4.0)
        A2(1, iRadius+1)=rx*rx/2.0+rx*dx2/3.0+dx2*dx2/12.0
        dx1=dx2
    END DO L_R

    dx1=rSpacings(1)
    dx2=rSpacings(nRSpacings)

    ! The special cases
    B(1, 1) = 0.0; B(2, 1)  = dx1*(r0+dx1/4.0)/3.
    A1(1, 1)= 0.0; A1(2, 1) = (r0+dx1/2.0)/dx1
    A2(1, 1)= 0.0; A2(2, 1) = -1.0*(r0*r0/2.0+r0*dx1/3.+dx1*dx1/12.0)

    B(2, nRadii)  = dx2*(r1-dx2/4.0)/3.0
    B(3, nRadii)  = 0.0
    A1(2, nRadii) = (r1-dx2/2.0)/dx2
    A1(3, nRadii) = 0.0
    A2(2, nRadii) = r1*r1/2.0-r1*dx2/3.0+dx2*dx2/12.0
    A2(3, nRadii) = 0.0

    RETURN
    END

    SUBROUTINE TUVWMatrices(iRun)
    USE ParallelTestData

    USE, INTRINSIC :: IEEE_Arithmetic

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: iRun
    real*8, DIMENSION(:, :), ALLOCATABLE :: Um, Vm, Wm

    real*8 :: r0, r1, dr, value
    INTEGER :: allocationStatus, deallocationStatus
    INTEGER :: i, k
    real*8, PARAMETER :: oneSixth=1.0/6, oneThird=1.0/3

    r0 = 5.9; r1=7.2
    dr  =  (r1-r0)/(nRadii-1)

    IF(ALLOCATED(Um)) THEN
        DEALLOCATE(Um, STAT = deallocationStatus)
    END IF
    ALLOCATE(Um(nRadii, 0:nRadii+1), STAT = allocationStatus)

    IF(ALLOCATED(Vm)) THEN
        DEALLOCATE(Vm, STAT = deallocationStatus)
    END IF
    ALLOCATE(Vm(nRadii, 0:nRadii+1), STAT = allocationStatus)

    IF(ALLOCATED(Wm)) THEN
        DEALLOCATE(Wm, STAT = deallocationStatus)
    END IF
    ALLOCATE(Wm(nRadii, 0:nRadii+1), STAT = allocationStatus)

    Um=0; Vm=0; Wm=0
    Vm(1,1) = r0/(2*dr) + oneSixth
    DO i=2,nRadii-1
        Vm(i,i) = (r0+(i-1)*dr)/dr
    END DO
    Vm(nRadii,nRadii) = r1/(2*dr) - oneSixth
    DO i=1,nRadii-1
        value = (r0+(i-1)*dr)/(2*dr) + oneSixth
        Um(i+1,i) = value
        Um(i,i) = -value
        Vm(i+1,i) = -value
        value = (r0+(i-1)*dr)/(2*dr) + oneThird
        Wm(i,i) = value
        Wm(i+1,i) = -value
        Vm(i,i+1) = -value
    END DO

    U=0; V=0; W=0
    DO i = 1, nRadii
        U(i, 1) = Um(i, i-1);  V(i, 1) = Vm(i, i-1);  W(i, 1) = Wm(i, i-1)
        U(i, 2) = Um(i, i);    V(i, 2) = Vm(i, i);    W(i, 2) = Wm(i, i)
        U(i, 3) = Um(i, i+1);  V(i, 3) = Vm(i, i+1);  W(i, 3) = Wm(i, i+1)
    END DO

    DEALLOCATE(Um, STAT = deallocationStatus)
    DEALLOCATE(Vm, STAT = deallocationStatus)
    DEALLOCATE(Wm, STAT = deallocationStatus)
    RETURN

    END SUBROUTINE TUVWMatrices

    subroutine delta_step (times)
    real*8         :: times(2)
!
    real*8         :: startTime
    integer*8      :: Startcount, count_rate
    real*8,save    :: lastTime  = 0
    integer*8,save :: lastcount = 0
!
      CALL SYSTEM_CLOCK (Startcount,count_rate)
      CALL CPU_TIME     (startTime)
!
      times(1) = times(1) + dble(Startcount-lastcount) / dble(count_rate)
      times(2) = times(2) + (startTime - lastTime)
      lastcount = Startcount
      lastTime  = startTime
    end subroutine delta_step
    