|
|
lurese (Programmer) |
3 Jun 12 16:27 |
Hi,
I use the ISML function DQDAGS to calculate integral. The function F1, F2 need to return arrays. I don't know how to deal with it. The following is part of the code.
thanks for every help.
lur
****************************
..............
DO 55 I=1, N_ZONE
A(I) = 0.078D0 + 0.02D0 * (I-1)
B(I) = 0.078D0 + 0.02D0 * (I)
55 ENDDO
DO 60 N = 1, N_TEXP
INTEGR1(N)=0.0D0
INTEGR2(N)=0.0D0
INTEGR3(N)=0.0D0
DO 70 I =1, N_ZONE
CALL DQDAGS (F1, A(I), B(I), ERRABS,ERRREL,ANS1(I, N),ERR1(I, N))
CALL DQDAGS (F2, A(I), B(I), ERRABS,ERRREL,ANS2(I, N),ERR2(I, N))
CALL DQDAGS (F3, A(I), B(I), ERRABS,ERRREL,ANS3(I, N),ERR3(I, N))
70 ENDDO
INTEGR1(N)= INTEGR1(N) + ANS1(I, N)
INTEGR2(N)= INTEGR2(N) + ANS2(I, N)
INTEGR3(N)= INTEGR3(N) + ANS3(I, N)
60 ENDDO
..............
C ****************************************************************
DOUBLE PRECISION FUNCTION F1(R, SPAREA, N_TEXP, N_ZONE)
DOUBLE PRECISION I, N, II, NN
DOUBLE PRECISION N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP)
DO 80 N = 1, N_TEXP
DO 81 I = 1, N_ZONE
F1(I, N) = SPAREA(I, N) * R
81 ENDDO
80 ENDDO
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
DOUBLE PRECISION FUNCTION F2(R, N_TEXP,N_ZONE, HOLLIQ, THFILM, SPAREA, DIADRO)
DOUBLE PRECISION I, N, II, NN
DOUBLE PRECISION N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION HOLLIQ(N_TEXP)
DOUBLE PRECISION THFILM(N_TEXP)
DOUBLE PRECISION DIADRO(N_ZONE, N_TEXP)
DO 90 N = 1, N_TEXP
DO 100 I =1, N_ZONE
F2(I, N) = 6.D0 * (HOLLIQ(N)- THFILM(N)) * R / DIADRO(I, N)
100 ENDDO
90 ENDDO
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
DOUBLE PRECISION FUNCTION F3(R)
DOUBLE PRECISION R !integration variable
F3= R
RETURN
END FUNCTION
C **************************************************************** |
|
You see how you are using a variable inside your function that has the same name as the function? When what you want to return is not just a scalar but an array, you need to actually, explicitly, declare such variable and specify it as an array with dimensions, etc. |
|
|
lurese (Programmer) |
3 Jun 12 16:52 |
Hi Salgerman,
I modify the code as following, but it still doesn't work.
****************************
..............
DO 55 I=1, N_ZONE
A(I) = 0.078D0 + 0.02D0 * (I-1)
B(I) = 0.078D0 + 0.02D0 * (I)
55 ENDDO
DO 60 N = 1, N_TEXP
INTEGR1(N)=0.0D0
INTEGR2(N)=0.0D0
INTEGR3(N)=0.0D0
DO 70 I =1, N_ZONE
CALL DQDAGS (F1, A(I), B(I), ERRABS,ERRREL,ANS1(I, N),ERR1(I, N))
CALL DQDAGS (F2, A(I), B(I), ERRABS,ERRREL,ANS2(I, N),ERR2(I, N))
CALL DQDAGS (F3, A(I), B(I), ERRABS,ERRREL,ANS3(I, N),ERR3(I, N))
70 ENDDO
INTEGR1(N)= INTEGR1(N) + ANS1(I, N)
INTEGR2(N)= INTEGR2(N) + ANS2(I, N)
INTEGR3(N)= INTEGR3(N) + ANS3(I, N)
60 ENDDO
..............
C ****************************************************************
DOUBLE PRECISION FUNCTION F1(R, SPAREA, N_TEXP, N_ZONE)
DOUBLE PRECISION I, N, II, NN
DOUBLE PRECISION N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP)
F1 = SPAREA(I, N) * R
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
DOUBLE PRECISION FUNCTION F2(R, N_TEXP,N_ZONE, HOLLIQ, THFILM, SPAREA, DIADRO)
DOUBLE PRECISION I, N, II, NN
DOUBLE PRECISION N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION HOLLIQ(N_TEXP)
DOUBLE PRECISION THFILM(N_TEXP)
DOUBLE PRECISION DIADRO(N_ZONE, N_TEXP)
F2 = 6.D0 * (HOLLIQ(N)- THFILM(N)) * R / DIADRO(I, N)
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
DOUBLE PRECISION FUNCTION F3(R)
DOUBLE PRECISION R !integration variable
F3= R
RETURN
END FUNCTION
C **************************************************************** |
|
Would be not too bad an idea, if you could be more specific on what did not work. Did you receive compiler-errors, runtime errors or is it you just did not receive the expected result ? Personally, I i am doubting if you can have functions to return arrays as a result. I think using subroutines would be better in this case. For instance, instead of your function f1 you could write: CODEsubroutine EvalF1 (F1, R, SPARAREA, N_TEXP, N_ZONE)
DOUBLE PRECISION F1 (N_TEXP, N_ZONE) ! declaration of f1 as array
DOUBLE PRECISION I, N, II, NN ! <----- you use these as integers !?!?
DOUBLE PRECISION N_TEXP, N_ZONE ! <----- you use these as integers !?!?
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP)
DO 80 N = 1, N_TEXP
DO 81 I = 1, N_ZONE
F1(I, N) = SPAREA(I, N) * R
81 ENDDO
80 ENDDO
RETURN
END You should rethink the types of your variables. Double precision as indices is not the best choice. Norbert The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true. |
|
lurese: I don't see where you declared your "function" variables to be arrays...for example, your original F1 function should look something like this: CODEC ****************************************************************
DOUBLE PRECISION FUNCTION F1(R, SPAREA, N_TEXP, N_ZONE)
double precision F1(N_ZONE, N_TEXP)
DOUBLE PRECISION I, N, II, NN
DOUBLE PRECISION N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP)
DO 80 N = 1, N_TEXP
DO 81 I = 1, N_ZONE
F1(I, N) = SPAREA(I, N) * R
81 ENDDO
80 ENDDO
RETURN
END FUNCTION |
|
Like Norbert said, you need to be more specific as to what exactly is not working...post the entire code if you need to, use "code" and "/code" tags enclosed each in square brackets (like these ones "[" and "]") or include it as an attachment if it is too large.
|
|
|
lurese (Programmer) |
4 Jun 12 11:25 |
According to the suggestions you gave, I post all the code here. Since the ISML function QDAGA require the function F1, F2, F3 to be the eternal, I can not declare it as a array. I modify the variable type with I, N, II, NN according to gummibaer's suggestion. Thanks. CODE --> EnglishC
PROGRAM INTERAREA
C ****************************************************************
C This is a code for calculating the effective interfacial area(ae)
C with multi zones in a reactor. Other detail information
C was given by the paper.
C luoyong980@gmail.com
C ****************************************************************
C ----------------------------------------------------------------
C Chapter 1. Variables declaration
C ----------------------------------------------------------------
IMPLICIT NONE
C *** Set number of packings ***
INTEGER N_PACK
PARAMETER (N_PACK = 8)
C *** Set number of experiments of each packing ***
INTEGER N_EXPR
PARAMETER (N_EXPR = 10)
C *** Set number of total experiments of all the packings ***
INTEGER N_TEXP
PARAMETER (N_TEXP = N_PACK * N_EXPR)
C *** Set number of packing zones in the rotor ***
INTEGER N_ZONE
PARAMETER (N_ZONE = 4)
C *** Set variables ***
INTEGER I,II
INTEGER N,NN
COMMON /DAT/ II, NN
COMMON /CONSTANTS/ DIAFIB, DENLIQ, EFAREA, FLOLIQ, RSPEED,
* SPAREA, TENLIQ, VISLIQ, RADIUS
C Original Input constant variables
DOUBLE PRECISION DENLIQ(N_TEXP) !Density of liquid, [Kg/m3]
DOUBLE PRECISION EFAREA(N_TEXP) !Effective interfacial area of..,ae [1/m]
DOUBLE PRECISION FLOLIQ(N_TEXP) !Liquid flow rate,L [L/h]
DOUBLE PRECISION RSPEED(N_TEXP) !Rotational speed, N [1/s]
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP) !Specfic surface area ..,ap [1/m]
DOUBLE PRECISION TENLIQ(N_TEXP) !Surface tension of the liquid, [Kg/s2]
DOUBLE PRECISION VISLIQ(N_TEXP) !Kinematic viscosity of liquid,v [m2/s]
DOUBLE PRECISION RADIUS(N_ZONE) !Radius of the rotor,r [1/m]
DOUBLE PRECISION DIAFIB(N_ZONE, N_TEXP) !Diameter of fiber,d [mm]
C Calculating variables (1 th)
DOUBLE PRECISION ACCELE(N_TEXP) !centrifugal acceleration,ac [m/s2]
DOUBLE PRECISION FLULIQ(N_TEXP) !Liquid flux, Lf [m3/(m2 s)]
C Calculating variables (2 th)
DOUBLE PRECISION HOLLIQ(N_TEXP) !Liquid holdup [m3/m3]
DOUBLE PRECISION THFILM(N_TEXP) !Thickness of liquid film [m]
DOUBLE PRECISION DIADRO(N_ZONE, N_TEXP) !Diameter of liquid droplet[m]
DOUBLE PRECISION PERCEN(N_TEXP) !Percentage of the surface area [-]
C Other variables
DOUBLE PRECISION HEIGHT !Height of the rotor,H [m]
DOUBLE PRECISION R !Integral variables
DOUBLE PRECISION B1(N_ZONE)
DOUBLE PRECISION B2(N_ZONE)
DOUBLE PRECISION INTEGR1(N_ZONE, N_TEXP) !Integral part 1
DOUBLE PRECISION INTEGR2(N_ZONE, N_TEXP) !Integral part 2
DOUBLE PRECISION INTEGR3(N_ZONE, N_TEXP) !Integral part 3
DOUBLE PRECISION SINTEGR1(N_TEXP) !Integral part 1
DOUBLE PRECISION SINTEGR2(N_TEXP) !Integral part 2
DOUBLE PRECISION SINTEGR3(N_TEXP) !Integral part 3
C Output results
OPEN (UNIT = 777, FILE = 'RESULT.out', STATUS = 'UNKNOWN')
C ----------------------------------------------------------------
C Chapter 2. Read input data
C ----------------------------------------------------------------
WRITE(6, *) 'READ INPUT DATA ........................'
OPEN (UNIT = 50, FILE = 'DIAFIB.INP', STATUS = 'OLD')
READ (50, *) ((DIAFIB(I, N), I = 1, N_ZONE), N = 1, N_TEXP)
CLOSE (UNIT = 50)
OPEN (UNIT = 55, FILE = 'DENLIQ.INP', STATUS = 'OLD')
READ (55, *) (DENLIQ(N), N = 1, N_TEXP)
CLOSE (UNIT = 55)
OPEN (UNIT = 60, FILE = 'EFAREA.INP', STATUS = 'OLD')
READ (60, *) (EFAREA(N), N = 1, N_TEXP)
CLOSE (UNIT = 60)
OPEN (UNIT = 65, FILE = 'FLOLIQ.INP', STATUS = 'OLD')
READ (65, *) (FLOLIQ(N), N = 1, N_TEXP)
CLOSE (UNIT = 65)
OPEN (UNIT = 70, FILE = 'RSPEED.INP', STATUS = 'OLD')
READ (70, *) (RSPEED(N), N = 1, N_TEXP)
CLOSE (UNIT = 70)
OPEN (UNIT = 75, FILE = 'SPAREA.INP', STATUS = 'OLD')
READ (75, *) ((SPAREA(I, N), I = 1, N_ZONE), N = 1, N_TEXP)
CLOSE (UNIT = 75)
OPEN (UNIT = 80, FILE = 'TENLIQ.INP', STATUS = 'OLD')
READ (80, *) (TENLIQ(N), N = 1, N_TEXP)
CLOSE (UNIT = 80)
OPEN (UNIT = 85, FILE = 'VISLIQ.INP', STATUS = 'OLD')
READ (85, *) (VISLIQ(N), N = 1, N_TEXP)
CLOSE (UNIT = 85)
OPEN (UNIT = 90, FILE = 'RADIUS.INP', STATUS = 'OLD')
READ (90, *) (RADIUS(N), N = 1, N_ZONE)
CLOSE (UNIT = 90)
WRITE(6, *) '................................... done'
WRITE(6, *) ''
C ----------------------------------------------------------------
C Chapter 3. Call subroutines to calculate the intermediate various
C ----------------------------------------------------------------
C ----------------------------------------------------------------
C Calculate integral 1, 2, 3
C ----------------------------------------------------------------
B1(N_ZONE)=0
B2(N_ZONE)=0
II=1
NN=1
DO 55 I=1, N_ZONE
B1(I) = 0.078D0 + 0.02D0 * (I-1)
B2(I) = 0.078D0 + 0.02D0 * (I)
55 ENDDO
DO 60 N = 1, N_TEXP
DO 70 I =1, N_ZONE
CALL MASSTR(N_ZONE, N_TEXP, B1(I), B2(I), INTEGR1(I, N),
* INTEGR2(I, N), INTEGR3(I, N))
II=II+1
70 ENDDO
NN=NN+1
60 ENDDO
DO 61 N = 1, N_TEXP
SINTEGR1(N)=0.0D0
SINTEGR2(N)=0.0D0
SINTEGR3(N)=0.0D0
DO 71 I =1, N_ZONE
SINTEGR1(N) = SINTEGR1(N)+ INTEGR1(I, N)
SINTEGR2(N) = SINTEGR2(N)+ INTEGR2(I, N)
SINTEGR3(N) = SINTEGR3(N)+ INTEGR3(I, N)
C WRITE(*,*) INTEGR1(I,N), INTEGR2(I,N), INTEGR3(I,N)
71 ENDDO
C WRITE(*,*) SINTEGR1(N), SINTEGR2(N), SINTEGR3(N)
61 ENDDO
C ----------------------------------------------------------------
C 6. Calculate the percentage
C ----------------------------------------------------------------
DO 5 N = 1, N_TEXP
PERCEN(N)= (EFAREA(N) * SINTEGR3(N) - SINTEGR2(N))/ SINTEGR1(N)
* *100.D0
5 END DO
C ----------------------------------------------------------------
C 7. Output the data
C ----------------------------------------------------------------
WRITE(777,100)
100 FORMAT(' No.',3X,'d1(mm)',3X,'d2(mm)',3X, 'd3(mm)',3X, 'd4(mm)',3X,
* 'N (1/s) ',3X,'Percentage (%)',3X)
WRITE(777,*)''
DO 10 N = 1, N_TEXP
WRITE(777,200) N, DIAFIB(1, N_TEXP), DIAFIB(2, N_TEXP),
* DIAFIB(3, N_TEXP), DIAFIB(4, N_TEXP), RSPEED(N), PERCEN(N)
200 FORMAT(I2, 3X, E10.2, 3X, E10.2, 3X, E10.2, 3X, E10.2, 3X,
* E10.2, 3X, E10.2, 3X)
10 END DO
WRITE(*,*)"Success!"
STOP
END
C ----------------------------------------------------------------
C MAIN PROGRAM STOPS HERE!!
C ----------------------------------------------------------------
C ----------------------------------------------------------------
C THE FOLLOWING ARE THE SUBROUTINES.
C ----------------------------------------------------------------
C ****************************************************************
SUBROUTINE ACCRPB(N_TEXP, ACCELE, RSPEED, R)
C ****************************************************************
C Routine name: ACCRPB
C Routine number: 1
C Modified by: Yong Luo at 05/30/2012
C
C Function: Calculate the centrifugal acceleration in a RPB
C
C ****************************************************************
INTEGER N, N_TEXP
DOUBLE PRECISION ACCELE(N_TEXP)!centrifugal acceleration,ac [m/s2]
DOUBLE PRECISION RSPEED(N_TEXP)!Rotational speed, N [1/s]
DOUBLE PRECISION R !integration variable
DO 20 N = 1, N_TEXP
ACCELE(N) = ((3.1415926D0* RSPEED(N) / 30.D0)**2.D0) * R
20 ENDDO
RETURN
END
C ****************************************************************
C ****************************************************************
SUBROUTINE HOLDUP(N_TEXP, HOLLIQ, ACCELE, FLULIQ, VISLIQ)
C ****************************************************************
C Routine name: HOLDUP
C Routine number: 2
C Modified by: Yong Luo at 05/30/2012
C
C Function: Calculates the liquid holdup
C
C ****************************************************************
INTEGER N, N_TEXP
DOUBLE PRECISION FLULIQ(N_TEXP) !Liquid flux, Lf [m3/(m2 s)]
DOUBLE PRECISION HOLLIQ(N_TEXP) !Liquid holdup [m3/m3]
DOUBLE PRECISION ACCELE(N_TEXP) !Centrifugal acceleration,ac [m/s2]
DOUBLE PRECISION VISLIQ(N_TEXP) !Kinematic viscosity of liquid,v [m2/s]
DO 30 N = 1, N_TEXP
HOLLIQ(N) =(0.039D0*(ACCELE(N)/100.0D0)**(-0.5D0)
* *(FLULIQ(N) /0.01D0)**0.6D0*(VISLIQ(N)/0.000001D0)
* **(0.22D0))**(1.0D0/1.6D0)
30 ENDDO
RETURN
END
C ****************************************************************
C ****************************************************************
SUBROUTINE THICKN(N_TEXP, THFILM, VISLIQ, FLULIQ, ACCELE)
C ****************************************************************
C Routine name: THICKN
C Routine number: 3
C Modified by: Yong Luo at 05/30/2012
C
C Function: Calculates the thickness of liquid films
C
C ****************************************************************
INTEGER N, N_TEXP
DOUBLE PRECISION FLULIQ(N_TEXP) !Liquid flux, Lf [m3/(m2 s)]
DOUBLE PRECISION VISLIQ(N_TEXP) !Kinematic viscosity of liquid,v [m2/s]
DOUBLE PRECISION ACCELE(N_TEXP) !centrifugal acceleration,ac [m/s2]
DOUBLE PRECISION THFILM(N_TEXP) !Thickness of liquid film [m]
c DOUBLE PRECISION SPAREA(N_TEXP) !Specfic surface area ..,ap [1/m]
DO 40 N = 1, N_TEXP
C !!!Note: THFILM= Liquid film thickness * specific surface area
THFILM(N) =4.2D0*(10.**8.0)*VISLIQ(N)*FLULIQ(N)
* / ACCELE(N)
40 ENDDO
RETURN
END
C ****************************************************************
C ****************************************************************
SUBROUTINE DIAMET(N_TEXP, N_ZONE, DIADRO, TENLIQ, DENLIQ,
* ACCELE, FLOLIQ, DIAFIB)
C ****************************************************************
C Routine name: DIAMET
C Routine number: 4
C Modified by: Yong Luo at 05/30/2012
C
C Function: Calculates the diameter of the liquid droplet
C
C ****************************************************************
INTEGER I, N, N_TEXP, N_ZONE
DOUBLE PRECISION ACCELE(N_TEXP) !Centrifugal acceleration; [m/s2]
DOUBLE PRECISION DIADRO(N_ZONE, N_TEXP) !Diameter of droplets,d [mm]
DOUBLE PRECISION DIAFIB(N_ZONE, N_TEXP) !Diameter of fiber,d [mm]
DOUBLE PRECISION DENLIQ(N_TEXP) !Density of liquid, [Kg/m3]
DOUBLE PRECISION FLOLIQ(N_TEXP) !Liquid flow rate,L [L/h]
DOUBLE PRECISION TENLIQ(N_TEXP) !Surface tension of the liquid, [Kg/s2]
DO 51 N = 1, N_TEXP
DO 52 I = 1, N_ZONE
DIADRO(I, N) =1.5158D0 * (TENLIQ(N) / (DENLIQ(N) * ACCELE(N)))
* **0.0495 * (FLOLIQ(N)/ 100.D0)**0.0735
* * (DIAFIB(I, N)/ 1.D0)**0.2489
52 ENDDO
51 ENDDO
RETURN
END
C ****************************************************************
C ****************************************************************
SUBROUTINE MASSTR(N_ZONE, N_TEXP, A1, A2, INTEGRAL1,
* INTEGRAL2, INTEGRAL3)
C ****************************************************************
C Routine name: MASSTR
C Routine number: 5
C Modified by: Yong Luo at 05/30/2012
C
C Function: Calculate the effective interfacial area
C
C ****************************************************************
USE IMSL
IMPLICIT NONE
INTEGER I, N, N_TEXP, N_ZONE,II,NN
DOUBLE PRECISION, EXTERNAL::F1, F2, F3
DOUBLE PRECISION INTEGRAL1
DOUBLE PRECISION INTEGRAL2
DOUBLE PRECISION INTEGRAL3
DOUBLE PRECISION A1
DOUBLE PRECISION A2
DOUBLE PRECISION, PARAMETER :: ERRABS = 1.E-5
DOUBLE PRECISION, PARAMETER :: ERRREL = 1.E-5
DOUBLE PRECISION ANS1
DOUBLE PRECISION ANS2
DOUBLE PRECISION ANS3
DOUBLE PRECISION ERR1
DOUBLE PRECISION ERR2
DOUBLE PRECISION ERR3
CALL DQDAGS (F1, A1, A2, ERRABS,ERRREL,ANS1,ERR1)
CALL DQDAGS (F2, A1, A2, ERRABS,ERRREL,ANS2,ERR2)
CALL DQDAGS (F3, A1, A2, ERRABS,ERRREL,ANS3,ERR3)
INTEGRAL1= ANS1
INTEGRAL2= ANS2
INTEGRAL3= ANS3
RETURN
END
C ****************************************************************
C ****************************************************************
REAL FUNCTION F1(R)
C ****************************************************************
C Routine name: F1
C Routine number: 6
C
C Function: integrand 1
C
C ****************************************************************
C *** Set number of packings ***
INTEGER N_PACK
PARAMETER (N_PACK = 8)
C *** Set number of experiments of each packing ***
INTEGER N_EXPR
PARAMETER (N_EXPR = 10)
C *** Set number of total experiments of all the packings ***
INTEGER N_TEXP
PARAMETER (N_TEXP = N_PACK * N_EXPR)
C *** Set number of packing zones in the rotor ***
INTEGER N_ZONE
PARAMETER (N_ZONE = 4)
COMMON /DAT/ II, NN
INTEGER I, N, II, NN
C INTEGER N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION SPAREA(N_ZONE, N_TEXP)!Specfic surface area ..,ap [1/m]
C DO 80 N = 1, N_TEXP
C DO 81 I = 1, N_ZONE
F1 = SPAREA(II, NN) * R
C F1 = R
C81 ENDDO
C80 ENDDO
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
REAL FUNCTION F2(R)
C ****************************************************************
C Routine name: F2
C Routine number: 7
C
C Function: integrand 2
C
C ****************************************************************
C *** Set number of packings ***
INTEGER N_PACK
PARAMETER (N_PACK = 8)
C *** Set number of experiments of each packing ***
INTEGER N_EXPR
PARAMETER (N_EXPR = 10)
C *** Set number of total experiments of all the packings ***
INTEGER N_TEXP
PARAMETER (N_TEXP = N_PACK * N_EXPR)
C *** Set number of packing zones in the rotor ***
INTEGER N_ZONE
PARAMETER (N_ZONE = 4)
COMMON /DAT/ II, NN
INTEGER I, N, II, NN
C INTEGER N_TEXP, N_ZONE
DOUBLE PRECISION R !integration variable
DOUBLE PRECISION HEIGHT
DOUBLE PRECISION FLULIQ(N_TEXP) !Liquid flux, Lf [m3/(m2 s)]
DOUBLE PRECISION HOLLIQ(N_TEXP) !Liquid holdup [m3/m3]
DOUBLE PRECISION ACCELE(N_TEXP) !Centrifugal acceleration,ac [m/s2]
DOUBLE PRECISION VISLIQ(N_TEXP) !Kinematic viscosity of liquid,v [m2/s]
DOUBLE PRECISION TENLIQ(N_TEXP) !Kinematic viscosity of liquid,v [m2/s]
DOUBLE PRECISION RSPEED(N_TEXP)!Rotational speed, N [1/s]
DOUBLE PRECISION FLOLIQ(N_TEXP)
DOUBLE PRECISION THFILM(N_TEXP) !Thickness of liquid film [m]
DOUBLE PRECISION DIADRO(N_ZONE, N_TEXP)!Diameter of droplets,d [mm]
DOUBLE PRECISION DIAFIB(N_ZONE, N_TEXP) !Diameter of fiber,d [mm]
DOUBLE PRECISION DENLIQ(N_TEXP) !Density of liquid, [Kg/m3]
C ----------------------------------------------------------------
C 1.Subroutine ACCRPB
C Function: Calculates the centrifugal acceleration in a RPB
C ----------------------------------------------------------------
CALL ACCRPB(N_TEXP, ACCELE, RSPEED, R)
C ----------------------------------------------------------------
C 2.Subroutine HOLDUP
C Function: Calculates the liquid holdup
C ----------------------------------------------------------------
HEIGHT = 0.05D0
DO 11 N=1, N_TEXP
FLULIQ(N) = (FLOLIQ(N)*0.001D0/3600.D0)
* /(2.D0*3.1415926D0* R * HEIGHT)
11 ENDDO
CALL HOLDUP(N_TEXP, HOLLIQ, ACCELE, FLULIQ, VISLIQ)
C ----------------------------------------------------------------
C 3.Subroutine THICKN
C Function: Calculates the thickness of liquid films
C ----------------------------------------------------------------
CALL THICKN(N_TEXP, THFILM, VISLIQ, FLULIQ, ACCELE)
C ----------------------------------------------------------------
C 4.Subroutine DIAMET
C Function: Calculates the diameter of the liquid droplet
C ----------------------------------------------------------------
CALL DIAMET(N_TEXP, N_ZONE, DIADRO, TENLIQ, DENLIQ,
* ACCELE, FLOLIQ, DIAFIB)
C ----------------------------------------------------------------
C DO 90 N = 1, N_TEXP
C DO 100 I =1, N_ZONE
F2 = 6.D0 * (HOLLIQ(NN)- THFILM(NN)) * R / DIADRO(II, NN)
C100 ENDDO
C90 ENDDO
RETURN
END FUNCTION
C ****************************************************************
C ****************************************************************
REAL FUNCTION F3(R)
C ****************************************************************
C Routine name: F3
C Routine number: 8
C
C Function: integrand 3
C
C ****************************************************************
DOUBLE PRECISION R !integration variable
F3= R
RETURN
END FUNCTION
C **************************************************************** |
|
|
lurese (Programmer) |
4 Jun 12 11:29 |
compiling information: CODE --> --------------------Configuration: interfacial area - Win32 Debug--------------------
Compiling Fortran...
D:\Fortran\1_Interfacial area\interfacial area2012-06-01\ae.for
D:\Fortran\1_Interfacial area\interfacial area2012-06-01\ae.for(413) : Warning: Variable SPAREA is used before its value has been defined
F1 = SPAREA(II, NN) * R
------------^
ae.obj - 0 error(s), 1 warning(s) CODE --> READ INPUT DATA ........................
................................... done
*** WARNING ERROR 2 from DQDAGS. Roundoff error has been detected. The
*** requested tolerances, ERRABS = 1.000000000000000D-05 and ERRREL
*** = 1.000000000000000D-05, cannot be reached.
*** WARNING ERROR 2 from DQDAGS. Roundoff error has been detected. The
*** requested tolerances, ERRABS = 1.000000000000000D-05 and ERRREL
*** = 1.000000000000000D-05, cannot be reached.
*** WARNING ERROR 2 from DQDAGS. Roundoff error has been detected. The
*** requested tolerances, ERRABS = 1.000000000000000D-05 and ERRREL
*** = 1.000000000000000D-05, cannot be reached.
*** WARNING ERROR 2 from DQDAGS. Roundoff error has been detected. The
*** requested tolerances, ERRABS = 1.000000000000000D-05 and ERRREL
*** = 1.000000000000000D-05, cannot be reached.
forrtl: severe (161): Program Exception - array bounds exceeded
Image PC Routine Line Source
interfacial area. 0040310B F1 409 ae.for
interfacial area. 00404B94 Unknown Unknown Unknown
interfacial area. 00403E80 Unknown Unknown Unknown
interfacial area. 00403BD9 Unknown Unknown Unknown
interfacial area. 00403862 Unknown Unknown Unknown
interfacial area. 00403032 MASSTR 361 ae.for
interfacial area. 004021C8 INTERAREA 146 ae.for
interfacial area. 004511D9 Unknown Unknown Unknown
interfacial area. 004366C9 Unknown Unknown Unknown
kernel32.dll 7C817077 Unknown Unknown Unknown
Press any key to continue |
|
lurese: please refer to my "3 Jun 12 19:50" post ...you have done nothing about that. Please read it and acknowledge that you understand it. |
|
lurese:
by the way, that source code as posted, shows that within the F1 function you keep overwriting the value of F1 with a scalar value, over and over and over within the nested do loops....you are totally missing something, here...again, refer to my post and suggested F1 calculation |
|
|
 |
Talk To Other Members
Read More Threads Like This One