Use DQDAGS to calculate integral 1

lurese · Jun 3, 2012

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 ****************************************************************

salgerman · Jun 3, 2012

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 · Jun 3, 2012

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 ****************************************************************

gummibaer · Jun 3, 2012

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:

Code:

    subroutine 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.

salgerman · Jun 3, 2012

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:

Code:

C     ****************************************************************
      DOUBLE PRECISION FUNCTION F1(R, SPAREA, N_TEXP, N_ZONE)
      [b]double precision F1(N_ZONE, N_TEXP)[/b]
      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

salgerman · Jun 3, 2012

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 · Jun 4, 2012

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:

C
      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 · Jun 4, 2012

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

salgerman · Jun 4, 2012

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.

salgerman · Jun 4, 2012

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

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Use DQDAGS to calculate integral 1

lurese

Programmer

salgerman

Programmer

lurese

Programmer

gummibaer

Programmer

salgerman

Programmer

salgerman

Programmer

lurese

Programmer

lurese

Programmer

salgerman

Programmer

salgerman

Programmer

Similar threads

Part and Inventory Search

Sponsor