!module xxx 
!contains

SUBROUTINE ROCKPROP(LPORE_STR,LPORE_DEP,STRNFRK,SIGAFRK,PORO_STR,PORO_DEP,DIAMGR_DEP,DIAMPT_STR,DIAMPT_DEP,DCRIT,NUBAR,DEL,LPORE,DIAM_GR)
REAL(KIND=8),INTENT(OUT) :: LPORE_STR,LPORE_DEP,STRNFRK,SIGAFRK,PORO_STR,PORO_DEP,DIAMGR_DEP,DIAMPT_STR,DIAMPT_DEP
REAL(KIND=8),INTENT(IN) :: DCRIT,NUBAR,DEL,LPORE,DIAM_GR
REAL(KIND=8) :: ASMP,BSMP,CSMP,ANS1,ANS2
integer, parameter:: OPTION_DCRIT=1
REAL(KIND=8),external :: FUNC1,FUNC2
real(kind=8) :: POROIN
POROIN=0.2
     !-----------------------------------------------------------------------------------------------------------
     IF (OPTION_DCRIT == 1) THEN 
         ! ASSUME CRITICAL PORE THROAT DIAMETER IS TWICE THE PARTICLE SIZE (TODD 1990)  
         ASMP    = 0.0D0
         BSMP    =-DLOG(DCRIT/NUBAR)/DEL
         CALL SIMPSON(ASMP, BSMP, FUNC1, ANS1)
         STRNFRK = 0.5 - ANS1
         SIGAFRK = 1.0 - STRNFRK
         PORO_STR   = STRNFRK*POROIN
         PORO_DEP   = POROIN - PORO_STR
         !-------------------------------------------------------------------------------------------------------
         ! CALCULATE THE LENGTH OF STRAINING PORES
         CSMP    = -100.0D0
         CALL SIMPSON(-BSMP, CSMP, FUNC1, ANS1)
         CALL SIMPSON(-BSMP, CSMP, FUNC2, ANS2)
         DIAMPT_STR = ANS2/ANS1
         LPORE_STR  = (DIAMPT_STR*LPORE)/(DIAM_GR*0.125)
         !-------------------------------------------------------------------------------------------------------       
         ! CALCULATE THE LENGTH AND GRAIN DIAMETER OF INTERCEPTION PORES.
         CSMP    = 100.0D0
         CALL SIMPSON(-BSMP, CSMP, FUNC1, ANS1)
         CALL SIMPSON(-BSMP, CSMP, FUNC2, ANS2)
         DIAMPT_DEP = ANS2/ANS1
         DIAMGR_DEP = DIAMPT_DEP/0.125
         LPORE_DEP  = DIAMGR_DEP*LPORE/DIAM_GR
     ELSEIF (OPTION_DCRIT == 2) THEN
        ! USE CIVAN'S MODEL TO CALCULATE CRITICAL PORE THROAT DIAMETER AS A FUNCTION OF PARTICLE REYNOLD NUMBER 
        OPEN (99,FILE = "WARNING FILE.lis")
        WRITE(99,*) "WARNING: PLEASE GIVE CIVAN'S MODEL TO CALCULATE CRITICAL PORE THROAT DIAMETER FUNCTION!"
        CLOSE(99)
        STOP
     ELSE
        OPEN (99,FILE = "WARNING FILE.lis")
        WRITE(99,*) "WARNING: PLEASE GIVE OTHER DEFINED CRITICAL PORE THROAT DIAMETER FUNCTION!"
        CLOSE(99)
        STOP
     ENDIF 
END SUBROUTINE ROCKPROP

FUNCTION FUNC1(X)
!----------------------------------------------------------------------------------------------------------------     
IMPLICIT NONE
REAL(KIND=8) :: X,FUNC1
real(kind=8) :: PI
PI=3.14
     !-----------------------------------------------------------------------------------------------------------
     FUNC1   = DEXP(-X**2.0D0/2.0D0)/(2.0D0*PI)**0.5D0
END FUNCTION FUNC1

FUNCTION FUNC2(X)
!----------------------------------------------------------------------------------------------------------------     
IMPLICIT NONE
REAL(KIND=8)  :: WEIGHT, DIMAPT, X,FUNC2
    REAL(KIND=8):: DCRIT,NUBAR,DEL,LPORE,DIAM_GR
    real(kind=8) :: PI
PI=3.14
    DCRIT=10.
    NUBAR=1.
    DEL  =0.5
    LPORE =1.
    DIAM_GR=100.
     !-----------------------------------------------------------------------------------------------------------
     WEIGHT  = DEXP(-X**2.0D0/2.0D0)/(2.0D0*PI)**0.5D0
     DIMAPT  = DEXP(X*DEL)*NUBAR
     FUNC2   = WEIGHT*DIMAPT
    END FUNCTION FUNC2

SUBROUTINE SIMPSON(ASMP, BSMP, FUNC, SUM)
REAL (KIND=8), INTENT(IN)  :: ASMP, BSMP
REAL (KIND=8) :: FUNC
REAL (KIND=8), INTENT(OUT) :: SUM
REAL (KIND=8) :: X, H, TOL, ERRACT, BOUND, S, SBAR, ERR, ERREST, ERRBND
REAL (KIND=8) :: F0, F25, F50, F75, F100
INTEGER       :: N
EXTERNAL FUNC
     !-----------------------------------------------------------------------------------------------------------
     ! INITIALIZE VARIABLE CONSTANTS
       IOUT  = 3
       TOL   = 0.5D-10
       INTMAX= 2000
       WRITE(IOUT,299)
       WRITE(IOUT,300) TOL
299    FORMAT(/,1X,'NUMERICAL INTEGRATION BY SIMPSOS RULE ADAPTIVE QUADRATURE',/)
300    FORMAT(/,1X, 'ERROR TOLERANCE : ', F15.13,/)
     !-----------------------------------------------------------------------------------------------------------
     ! SELECT APPROPRIATE INITIAL NUMBER OF SUBINTERVALS DEPENDING ON ERROR TOLERANCE
       N    = 32
       H    = (BSMP-ASMP)/N
     !-----------------------------------------------------------------------------------------------------------
     ! INITIALIZE RUNNING SUMS AND INTEGER VARIABLES
       SUM    = 0.D0
       ERREST = 0.D0
       KCNT   = 0
       ITEX   = 0
       NOPROG = 0
       BOUND  = 0.0D0
       X      = ASMP
       F0     = FUNC(ASMP)
11     CONTINUE
       F25    = FUNC(X+H/4.0D0)
       F50    = FUNC(X+H/2.0D0)
       F75    = FUNC(X+3.0D0*H/4.0D0)
       F100   = FUNC(X+H)
12     CONTINUE
     !-----------------------------------------------------------------------------------------------------------
     ! CALCULATE ONE-PANEL AND TWO-PANEL SUMS
       S      = H/6.0D0*(F0 + 4.0D0*F50 + F100)
       SBAR   = H/12.0D0*(F0+4.0D0*F25+2.0D0*F50+4.0D0*F75+F100)
     !-----------------------------------------------------------------------------------------------------------
     ! TERMINATE INTEGRATION IF NO PROGRESS IS BEING MADE OR INTMAX IS EXCEEDED
       IF(SUM.EQ.SUM+SBAR) THEN
           NOPROG = 1
           CONTINUE
       ELSEIF(N.GT.INTMAX) THEN
           ITMAX  = 1
           CONTINUE
       ENDIF
     !-----------------------------------------------------------------------------------------------------------
     ! UPDATE THE SUMS IF TOL IS SATISFIED. IF NOT, TRY AGAIN
       ERR  = (SBAR-S)/15.0D0 
       IF (DABS(ERR) .LT. TOL*H/(BSMP-ASMP) .OR. NSD .GT. 5) THEN
           SUM = SUM+SBAR
           X     = X+H
           ERREST= ERREST+ERR
           IF(NSD .GT. 5) KCNT = KCNT + 1
           NSD   = 0
           IF (X .GE. BSMP) THEN
               X = X - H
               GO TO 100
           ENDIF
           N    = N + 1
           F0   = F100
           GO TO 11
       ELSE
           H    = H/2.0D0
           F100 = F50
           F50  = F25
           F25  = FUNC(X + H/4.0D0)
           F75  = FUNC(X+3.0D0*H/4.0D0)
           NSD  = NSD + 1
           GO TO 12
       ENDIF
100    CONTINUE
     !-----------------------------------------------------------------------------------------------------------
     ! PRINT THE RESULTS
       WRITE(IOUT, 302) ASMP, X, N, SUM
302    FORMAT(1X,'INTEGRAL FROM', E18.11, 'TO ', E18.11, ' OVER ', I4, ' POINTS= ', F21.18)
       WRITE(IOUT,303) ERREST/SUM
303    FORMAT(1X, 'ESTIMATED RELATIVE ERROR : ', E16.7)
       IF(NOPROG .EQ. 1) WRITE(IOUT,308)
       IF(ITEX .EQ. 1)   WRITE(IOUT,309)
308    FORMAT(1X,'INTEGRATION STOPPED - NO PROGRESS BEING MADE')
309    FORMAT(1X,'INTEGRATION STOPPED - MAXIMUM ITERATIONS EXCEEDED')
       IF(KCNT .GT. 0) WRITE(IOUT,310) KCNT
310    FORMAT(1X,'CONVERGENCE FAILED', I4, 'TIMES')
END SUBROUTINE SIMPSON 
    !end module xxx
!================================================================================================================
program testprog
    !use xxx
    REAL(KIND=8):: LPORE_STR,LPORE_DEP,STRNFRK,SIGAFRK,PORO_STR,PORO_DEP,DIAMGR_DEP,DIAMPT_STR,DIAMPT_DEP
    REAL(KIND=8):: DCRIT,NUBAR,DEL,LPORE,DIAM_GR
    DCRIT=10.
    NUBAR=1.
    DEL  =0.5
    LPORE =1.
    DIAM_GR=100.
    CALL ROCKPROP(LPORE_STR,LPORE_DEP,STRNFRK,SIGAFRK,PORO_STR,PORO_DEP,DIAMGR_DEP,DIAMPT_STR,DIAMPT_DEP,DCRIT,NUBAR,DEL,LPORE,DIAM_GR)
    write(*,*) LPORE_STR,LPORE_DEP,STRNFRK,SIGAFRK,PORO_STR,PORO_DEP,DIAMGR_DEP,DIAMPT_STR,DIAMPT_DEP
pause
end program testprog