Deep Earth Physics - Dipolarityを計算するプログラム

Dipolarityを計算するプログラム

この間まで載せていたのはすごく間違っていたので作り直した。

↑今まで、プログラムを組む前の式をノートに書いていて、検討のたびに書き直して色んなところに散逸してしまっていたので、ちゃんとまとめてみました。プログラムはこの通りに作ったつもりです。

program dipolarity
!! 12次までのガウス係数@CMB を読み込んで色々計算するプログラム !
implicit none
integer,parameter ::tmax=5000*17,tmin=1,nl=12 ! input patameters integer :: t,i,j,delta integer, dimension(:,:,:), allocatable ::l,m real(8), dimension(:,:,:,:), allocatable :: vout real(8), dimension(:,:,:), allocatable :: glm,hlm real(8), dimension(:), allocatable :: time,fdip1,fdip2,br2_sphere,br2_sphere1,br2ds,br2ds1 real(8), parameter :: r_core=330.0d3,r_moon=1737.0d3,eta=0.2d0 real(8), parameter :: pi=3.14159265358979323846264338327950288d0 real(8),parameter :: pm=5.0d0 real(8) :: nm,r_o
allocate(l(1:tmax,0:nl,0:nl),m(1:tmax,0:nl,0:nl)) allocate(vout(1:tmax,0:1,0:nl,0:nl)) allocate(time(1:tmax),fdip1(1:tmax),fdip2(1:tmax)) allocate(br2ds(1:tmax),br2ds1(1:tmax),br2_sphere(1:tmax),br2_sphere1(1:tmax)) allocate(glm(1:tmax,0:nl,0:nl),hlm(1:tmax,0:nl,0:nl))
r_o=1.0d0/(1.0d0-eta)
!nmの設定 do i=0,nl if(m(t,j,i)==0) then nm=2.0d0*pi else nm=pi endif enddo
!deltaの設定 do i=0,nl if(m(t,j,i)==0) then delta=1 else delta=0 endif enddo
!ファイルを読み込む open(1,file="fort.20",status="old") do t=1,tmax read(1,*) time(t) write(*,*) time(t) do i=0,1 do j=i,nl read(1,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(1)
open(2,file="fort.21",status="old") do t=1,tmax read(2,*) time(t) write(*,*) time(t) do i=2,3 do j=i,nl read(2,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(2)
open(3,file="fort.22",status="old") do t=1,tmax read(3,*) time(t) write(*,*) time(t) do i=4,5 do j=i,nl read(3,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(3)
open(4,file="fort.23",status="old") do t=1,tmax read(4,*) time(t) write(*,*) time(t) do i=6,7 do j=i,nl read(4,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(4)
open(7,file="fort.24",status="old") do t=1,tmax read(7,*) time(t) write(*,*) time(t) do i=8,9 do j=i,nl read(7,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(7)
open(8,file="fort.25",status="old") do t=1,tmax read(8,*) time(t) write(*,*) time(t) do i=10,11 do j=i,nl read(8,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(8) open(9,file="fort.26",status="old") do t=1,tmax read(9,*) time(t) write(*,*) time(t) do i=12,12 do j=i,nl read(9,*) l(t,j,i),m(t,j,i),vout(t,0,j,i),vout(t,1,j,i) enddo enddo enddo close(9)
!時間ごとにDipolarityの計算をする open(10,file="fdip.dat",status="replace") do t=1,tmax do i=0,nl do j=i,nl br2_sphere(t) = br2_sphere(t) && + ( dble(l(t,j,i))*dble(l(t,j,i))*(dble(l(t,j,i))+1.0d0)*(dble(l(t,j,i))+1.0d0) && * nm * ( vout(t,0,j,i)*vout(t,0,j,i) + vout(t,1,j,i)*vout(t,1,j,i) )/(r_o*r_o) ) if(l(t,j,i)==1) then br2_sphere1(t) = br2_sphere1(t) && + ( dble(l(t,j,i))*dble(l(t,j,i))*(dble(l(t,j,i))+1.0d0)*(dble(l(t,j,i))+1.0d0) && * nm * ( vout(t,0,j,i)*vout(t,0,j,i) + vout(t,1,j,i)*vout(t,1,j,i) )/(r_o*r_o) ) endif enddo enddo fdip1(t)=dsqrt(br2_sphere1(t)/br2_sphere(t)) write(10,*) time(t), fdip1(t) enddo close(10)

!ガウス係数からDipolarityを求める open(11,file="fdip_gauss.dat",status="replace")
do t=1,tmax do i=0,nl do j=i,nl glm(t,j,i) = dble(l(t,j,i))/(r_o*r_o) * ((r_core/r_moon)**(dble(l(t,j,i))+2.0d0)) && * dsqrt((2.0d0*dble(l(t,j,i))+1.0d0)/(2.0d0*(2.0d0-delta))) * (vout(t,0,j,i))
hlm(t,j,i) = dble(l(t,j,i))/(r_o*r_o) * ((r_core/r_moon)**(dble(l(t,j,i))+2.0d0)) && * dsqrt((2.0d0*dble(l(t,j,i))+1.0d0)/(2.0d0*(2.0d0-delta))) * (vout(t,1,j,i)) enddo enddo enddo
do t=1,tmax do i=0,nl do j=i,nl br2ds(t) = br2ds(t) && + ( ((dble(l(t,j,i))+1.0d0)*(dble(l(t,j,i))+1.0d0) ) && *( (r_moon/r_core)**(2.0d0*dble(l(t,j,i))+2.0d0) ) && *( glm(t,j,i)*glm(t,j,i) + hlm(t,j,i)*hlm(t,j,i) ) && *(2.0d0*(2.0d0-delta)/(2.0d0*dble(l(t,j,i))+1.0d0)) )*nm/(r_o*r_o)
if(l(t,j,i)==1) then br2ds1(t) = br2ds1(t) && + ( ((dble(l(t,j,i))+1.0d0)*(dble(l(t,j,i))+1.0d0) ) && *( (r_moon/r_core)**(2.0d0*dble(l(t,j,i))+2.0d0) ) && *( glm(t,j,i)*glm(t,j,i) + hlm(t,j,i)*hlm(t,j,i) ) && *(2.0d0*(2.0d0-delta)/(2.0d0*dble(l(t,j,i))+1.0d0)) )*nm/(r_o*r_o) endif enddo enddo fdip2(t)=dsqrt(br2ds1(t)/br2ds(t)) write(11,*) time(t), fdip2(t) enddo close(11)
end program dipolarity

よさそうな気がする．

あとは平均をとればよい．上の場合だと，だいたいfdip~0.45くらいになる気がする．