我编了一个利用IMSL库函数调用二重积分子函数TWODQ的程序，也就是说利用TWODQ程序分别对F，FIM进行积分，但是出了点问题。第一个函数F没问题，但是第二个被积函数FIM就不行了。请各路神仙帮忙看看好吗
    program main
    use IMSL
    INTEGER IRULE
    REAL A,B,ERRABS,ERRREL,RESULT,ERREST
    REAL G,H,F
    COMMON /PI/PI
    EXTERNAL G,H,F,FIM
    A=0.
    B=12.
    ERRABS=0.
    ERRREL=0.01
    IRULE=6
    CALL TWODQ(F,A,B,G,H,ERRABS,ERRREL,IRULE,RESULT,ERREST)
    WRITE(*,*)RESULT,ERREST

    CALL TWODQ(FIM,A,B,G,H,ERRABS,ERRREL,IRULE,RESULT,ERREST)
    
    WRITE(*,*)RESULT,ERREST
    END
***********************************
    BLOCK DATA
    REAL PI
    COMMON /PI/PI
    DATA PI/3.1415926/
    END BLOCK DATA
*********************************
    REAL FUNCTION F(R,THE)
    REAL R,THE,MASS
    COMMON /PI/PI

    Q=0.5

    RHO0=0.153
    A0=0.8
    BETA2=-0.113
    BETA4=-0.027
    Y20=0.25*SQRT(5./PI)*(-1+3*COS(THE)**2)
    Y40=(3./(16.*SQRT(PI)))*(3-30*COS(THE)**2+35*COS(THE)**4)
    MASS=131.
    R0=1.1
    RTHE=R0*(MASS**(1./3.))*(1+BETA2*Y20+BETA4*Y40)
    RHO=RHO0/(1+EXP((R-RTHE)/A0))

    F=2*PI*RHO*COS(Q*R*COS(THE))*(R**2)*SIN(THE)/MASS

    END
*************************************************************
    REAL FUNCTION FIM(R,THE)
    REAL R,THE,MASS
    COMMON /PI/PI

    Q=0.5

    RHO0=0.153
    A0=0.8
    BETA2=-0.113
    BETA4=-0.027
    Y20=0.25*SQRT(5./PI)*(-1+3*COS(THE)**2)
    Y40=(3./(16.*SQRT(PI)))*(3-30*COS(THE)**2+35*COS(THE)**4)
    MASS=131.
    R0=1.1
    RTHE=R0*(MASS**(1./3.))*(1+BETA2*Y20+BETA4*Y40)
    RHO=RHO0/(1+EXP((R-RTHE)/A0))

    FIM=2*PI*RHO*SIN(Q*R*COS(THE))*(R**2)*SIN(THE)/MASS
    END

    REAL FUNCTION H(X)
    COMMON/PI/PI
    H=PI
    END    

    REAL FUNCTION G(X)
    G=0.
    END