!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