[求助] 编译通过, 一直运行中, 函数调用无法返回!

program test
    
    implicit none
    
    integer :: i=0

    real(kind=4),external :: fun_real_g
    
    real(kind=4),dimension(10) :: x=(/10.0, 100.0, 200.0, 300.0, &
    400.0, 500.0, 600.0, 700.0,800.0, 900.0 /) 

    do i=1,10,1
        
        write(*,*) 'in test!!!!!!'
        
     write(*,*) i, x(i), fun_real_g( x(i) )
    
    end do
    
end program test 



  

real(kind=4) function fun_real_g ( r )
    

    implicit none
    

    
    real(kind=4),intent(in) :: r                       !计算点处的半径
    
    real(kind=4),external :: mass_of_sphere

    integer,parameter :: rho_layer_number=9, poly_coef=4 
    
    !密度分层数(rho layer number)为9, 每层均用3次(order)多项式(polynomial)来表达, 系数共有4个.
    
    real, parameter :: PI=3.141593
    real(kind=4), parameter :: there_epsilon=10.0E-30
    
    integer(kind=4) :: which_layer=0    !层号判断    
    integer(kind=4) :: i=0               !循环标    
    
    
    real(kind=4), dimension(poly_coef,rho_layer_number) :: rho_matrix
    real(kind=4), dimension(rho_layer_number) :: radius_matrix
    
!    用4行9列的矩阵来表达密度

    rho_matrix(:,1)= (/13.0885,  0.0   , -8.8381,  0.0   /)
    rho_matrix(:,2)=(/12.5815, -1.2638, -3.6426, -5.5281/)
    rho_matrix(:,3)=(/ 7.9565, -6.4761, +5.5283, -3.0807/)
    rho_matrix(:,4)=(/ 5.3197, -1.4836,  0.0   ,  0.0   /)
    rho_matrix(:,5)=(/11.2494, -8.0298,  0.0   ,  0.0   /)
    rho_matrix(:,6)=(/ 7.1089, -3.8045,  0.0   ,  0.0   /)
    rho_matrix(:,7)=(/ 2.6910, +0.6924,  0.0   ,  0.0   /)
    rho_matrix(:,8)=(/ 2.900,   0.0,     0.0   ,  0.0   /)
    rho_matrix(:,9)=(/ 2.600,   0.0,     0.0   ,  0.0   /)
    
    
    
!    用9行1列的矩阵来表达密度间断面的位置

    radius_matrix=(/1221.5, 3480.0, 5701.0,&
                    5771.0, 5971.0, 6151.0,&
                    6346.6, 6356.0, 6371.0/)
                    
    fun_real_g = 0.0 !
    
    write(*,*) rho_matrix
    write(*,*) radius_matrix
    
    
    
    write(*,*) 'in fun_real_g',fun_real_g

    
    if (r<=there_epsilon) then
        
        write(*,*) 'in fun_real_g line 66 !!'
        
        which_layer=0
        
        fun_real_g = 0.0
        
    else if (r<=1221.5) then
        which_layer=1
!         Inner core        
    else if (r<=3480.0) then
        which_layer=2
!         Outer core        
    else if (r<=5701.0) then
        which_layer=3
!         Lower mantle part 1, 2 and 3        
    else if (r<=5771.0) then
        which_layer=4
!         Transition zone part1
    else if (r<=5971.0) then
        which_layer=5
!         Transition zone part2        
    else if (r<=6151.0) then
        which_layer=6
!         Transition zone part3        
    else if (r<=6436.6) then
        which_layer=7
!          LVZ (low velocity zone) and LID (lid of LVZ)        
    else if (r<=6356.0) then
        which_layer=8
!         Lower crust        
    else if (r<=6371.0) then
        which_layer=9
!         Upper crust (本来还有一个海水层, 在此忽略!)        
    else    
        write(*,*) "radius given by user is wrong!"
        write(*,*) "radius of Earth Model should not be larger than 6371.0km!!!"
    end if

    
 
    do while(which_layer>=1)

        if (which_layer>=2) then  
        
            do i = 1, which_layer-1, 1
                fun_real_g = mass_of_sphere( (rho_matrix(:,i)-rho_matrix(:,i+1)) , radius_matrix(i) )
            end do
        
        end if
        
        fun_real_g = fun_real_g +  mass_of_sphere( rho_matrix(:,which_layer) , r )      

    end do  


end function fun_real_g


real(kind=4) function mass_of_sphere(rho_coeffs, r)
    
!        当rho_coeffs是密度多项式的系数    A0,A1,...时, rh(r/R) = A0 + A1*(r/R) +  A2*(r/R)^2 + ... , 其中R是规化半径;
!        一个球体的质量表达为: 4*PI*r^3* (  A0/3 + A1/4*(r/R) +  A2/5*(r/R)^2 + A3/6*(r/R)^3 + ...  ).
!        从而, 其质量系数为: A0/3, A1/4, A2/5, ...
!        调用horner来计算多项式!
        
    implicit none
    real(kind=4), dimension(4), intent(in) :: rho_coeffs
    real(kind=4), intent(in) :: r
    real(kind=4),external :: horner

    real(kind=4), parameter :: R_earth=6371.0
    real(kind=4), parameter :: PI=3.14159265
    
    real(kind=4),dimension(4) :: mass_ceffs
    
    integer(kind=4) :: j=0
    
    do j = 1, 4, 1    
        mass_ceffs(j) = rho_coeffs(j) / ( real( ( j+3),kind=4 ) )    
    end do
    
    mass_of_sphere= 4.0*PI*r**3*( horner(mass_ceffs, r/R_earth) )

end function mass_of_sphere
    
    

real(kind=4)  function horner (coeffs, x)
 
    implicit none
 
    real(kind=4), dimension (4), intent (in) :: coeffs
    real(kind=4), intent (in) :: x

    integer(kind=4) :: i
    
    horner = 0.0
    
    do i = 4, 1, -1
        horner = horner * x + coeffs (i)
    end do   
    

end function horner