RungeKutta.f90

Go to the documentation of this file.

MODULE rungekutta
  
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE rk4sum(numDim,numPoints,bufferSize,bufferInterval,h,K1,K2,K3,K4,stateData)
    
    INTEGER(KIND=4), INTENT(IN)      :: numDim
    INTEGER(KIND=8), INTENT(IN)      :: numPoints,bufferSize(numdim)
    INTEGER(KIND=8), INTENT(IN)      :: bufferInterval(2*numdim)
    REAL(KIND=8),    INTENT(IN)      :: h
    REAL(KIND=8),    INTENT(IN)      :: K1(numpoints)
    REAL(KIND=8),    INTENT(IN)      :: k2(numpoints)
    REAL(KIND=8),    INTENT(IN)      :: K3(numpoints)
    REAL(KIND=8),    INTENT(IN)      :: k4(numpoints)
    REAL(KIND=8),    INTENT(INOUT)   :: stateData(numpoints)
    
    INTEGER(KIND=8) :: I, J, K
    INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
    INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
    REAL(KIND=8)    :: fac1,fac2
  
    fac1 = h/6.0_8
    fac2 = h/3.0_8
    
    istart = bufferinterval(1)
    iend = bufferinterval(2)
    
    xsize = buffersize(1)
    
    IF(numdim == 3) THEN
       
       jstart = bufferinterval(3)
       jend   = bufferinterval(4)
       kstart = bufferinterval(5)
       kend   = bufferinterval(6)
       
       nplane = xsize*buffersize(2)
       
       DO k = kstart, kend
          zindex = (k-1)*nplane
          DO j = jstart, jend
             yzindex = (j-1)*xsize + zindex
             DO i = istart, iend
                bufferindex = yzindex + i
                statedata(bufferindex) = fac1*(k1(bufferindex)+k4(bufferindex)) + &
                     fac2*(k2(bufferindex) + k3(bufferindex)) + statedata(bufferindex)
             END DO
          END DO
       END DO
       
    ELSE IF (numdim == 2) THEN
       
       jstart = bufferinterval(3)
       jend   = bufferinterval(4)
       
       DO j = jstart, jend
          yindex = (j-1)*xsize
          DO i = istart, iend
             bufferindex = yindex + i
             statedata(bufferindex) = fac1*(k1(bufferindex)+k4(bufferindex)) + &
                  fac2*(k2(bufferindex) + k3(bufferindex)) + statedata(bufferindex)
          END DO
       END DO
       
    ELSE IF (numdim == 1) THEN
       
       DO bufferindex = istart, iend
          statedata(bufferindex) = fac1*(k1(bufferindex)+k4(bufferindex)) + &
               fac2*(k2(bufferindex) + k3(bufferindex)) + statedata(bufferindex)
       END DO
       
    ENDIF
    
    
  END SUBROUTINE rk4sum
  
END MODULE rungekutta