[Fortran] 纯文本查看 复制代码
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