!  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(10000,3) :: A1, B, A2
    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
    INTEGER :: iRun, iRadius, iSpecies, nThreads, nRSpacings, iDt
    REAL(8), EXTERNAL :: DiffCoeff, SedCoeff, ZCalc
    REAL(8), DIMENSION(0:10000+1) :: Dradius,sRadius
    REAL(8), DIMENSION(10000) :: Z
    REAL(8), DIMENSION(10000+1,3) :: concTotal
    REAL(8), DIMENSION(10000,3) :: G
    REAL(8), DIMENSION(3) :: rSpacings
    REAL(8) :: startTime, endTime, CPUTime, elsapsedTime, sum
    real(8) :: times(2,5)
    INTEGER :: countStart,count_rate,countEnd, op_n
    CHARACTER(LEN=100) :: msg

    integer(8) :: count
!
    nThreads=-1
!
    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.

    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))
      times = 0
      
 CALL SYSTEM_CLOCK (countStart,count_rate)
 CALL CPU_TIME     (startTime)
    sum   = 0
    count = 0
    DO iDt=1,100000
        dt=dt*1.00002
!
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP  DO PRIVATE(iRadius) schedule(static)
        DO iRadius = 1, nRadii
            Dradius(iRadius) = DiffCoeff (iSpecies, iRadius, concTotal(iRadius, 1:3)) ! diffusion coefficient at iRadius
            sRadius(iRadius) = SedCoeff  (iSpecies, iRadius, concTotal(iRadius, 1:3)) ! if compression, this will account for it
        END DO
!$OMP  END DO
!$OMP END PARALLEL

      call delta_step (times(1,1))

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP  DO PRIVATE(iRadius) schedule(static)

        L_R2_Parallel: &
            DO iRadius = 1, nRadii
                Z(iRadius)   = ZCalc(iRadius, iSpecies)
                G(iRadius, 1)= Dradius(iRadius-1)*dt*A1(iRadius, 1)+B(iRadius, 1)-sRadius(iRadius-1)*omSqRun*dt*A2(iRadius, 1)
                G(iRadius, 2)= Dradius(iRadius)  *dt*A1(iRadius, 2)+B(iRadius, 2)-sRadius(iRadius)  *omSqRun*dt*A2(iRadius, 2)
                G(iRadius, 3)= Dradius(iRadius+1)*dt*A1(iRadius, 3)+B(iRadius, 3)-sRadius(iRadius+1)*omSqRun*dt*A2(iRadius, 3)
            END DO L_R2_Parallel
!$OMP  END DO
!$        nThreads=omp_get_num_threads()
!$OMP END PARALLEL
!
      call delta_step (times(1,2))
      sum = sum                                               &
          + DOT_PRODUCT(Z(1:nRadii),   G(1:nRadii,1))         &
          + DOT_PRODUCT(G(1:nRadii,2), G(1:nRadii,3))

      call delta_step (times(1,3))
      op_n = 0

!$OMP PARALLEL DO                     &
!$OMP   DEFAULT  (SHARED)             &
!$OMP   PRIVATE  (iRadius)            & 
!$OMP   schedule (static)             & 
!$OMP   REDUCTION(+ : op_n)

            DO iRadius = 1, nRadii
              op_n = op_n+1
            END DO

!$OMP END PARALLEL DO
!
      call delta_step (times(1,4))
      count = count+1

      call delta_step (times(1,5))
    END DO
    
    CALL SYSTEM_CLOCK(countEnd)
    CALL CPU_TIME    (endTime)
    
    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
    WRITE(*,"(A)") msg
    !READ(*,*)
    
    write (*,*) ' '
    write (*,*) '!$omp entries',count
    write (*,*) 'omp loop size',op_n
    write (*,11) 'loop1', times(:,1), times(1,1)/dble(count), '  function calls'
    write (*,11) 'loop2', times(:,2), times(1,2)/dble(count), '  calculation'
    write (*,11) 'sum  ', times(:,3), times(1,3)/dble(count), '  dot_product'
    write (*,11) 'loop3', times(:,4), times(1,4)/dble(count), '  OMP overhead'
    write (*,11) 'timer', times(:,5), times(1,5)/dble(count), '  timer overhead'
 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(iRadius, 1)*Cm(iRadius-1, iSpecies)   &
        +B(iRadius, 2)*Cm(iRadius, iSpecies)     &
        +B(iRadius, 3)*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(iRadius, 1)*Ca(iRadius-1) & ! this part is for s
        + A2(iRadius, 2)*Ca(iRadius)   &
            + A2(iRadius, 3)*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(10000+1, 3), 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(iRadius, 2)=(dx1+dx2)*(rx/3.-(dx1-dx2)/12.)
        B(iRadius, 3)=dx2*(rx+dx2/2.0)/6.0
        B(iRadius+1, 1)=B(iRadius, 3)

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

        A2(iRadius, 2)=((dx1*dx1-dx2*dx2)/4. - rx*(dx1+dx2))/3.0
        A2(iRadius, 3)=-1.0*(rx*rx/2.0+2.0*rx*dx2/3.0+dx2*dx2/4.0)
        A2(iRadius+1, 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(1, 2)  = dx1*(r0+dx1/4.0)/3.
    A1(1, 1)= 0.0; A1(1, 2) = (r0+dx1/2.0)/dx1
    A2(1, 1)= 0.0; A2(1, 2) = -1.0*(r0*r0/2.0+r0*dx1/3.+dx1*dx1/12.0);

    B(nRadii, 2)  = dx2*(r1-dx2/4.0)/3.0; B(nRadii, 3)  = 0.0
    A1(nRadii, 2) = (r1-dx2/2.0)/dx2; A1(nRadii, 3) = 0.0
    A2(nRadii, 2) = r1*r1/2.0-r1*dx2/3.0+dx2*dx2/12.0; A2(nRadii, 3) = 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
    