program egy! implicit none
character(len=2) :: var1="bp" character(len=2) :: var2="bt" character(len=1) :: var="b" real(8),parameter :: eta=0.5d0
! real(8),parameter :: ek=1.2347172440720542E-004 !xi=0.1 ! real(8),parameter :: ek=1.5625002604693504E-004 !xi=0.2 ! real(8),parameter :: ek=2.0408282070913058E-004 !xi=0.3 ! real(8),parameter :: ek=2.7777970975710912E-004 !xi=0.4 real(8),parameter :: ek=4.0000027713744306E-004 !xi=0.5 ! real(8),parameter :: ek=6.2500044844826378E-004 !xi=0.6 ! real(8),parameter :: ek=1.1111114976624307E-003 !xi=0.7
real(8), parameter :: r_core=330.0d3,r_moon=1737.0d3 real(8),parameter :: pm=5.0d0 real(8),parameter :: pi=3.14159265358979323846264338327950288d0 real(8) :: time,ro,ri,volume,nm,egy_r,egy_t,egy_p,ramda,egy_cmb,egy_surface real(8) :: b_r2ds,b_theta2ds,b_phi2ds,oct,qua,dip real(8) :: b_r2ds1,b_theta2ds1,b_phi2ds1,b2ds,b2ds1 integer,parameter :: nr=99,nl=127,nt=192,np=400 integer :: k,kk,mm,ll,delta
integer,dimension(:,:),allocatable :: l,m integer,dimension(:),allocatable :: mp,lp real(8),dimension(:),allocatable :: m_egy_r,m_egy_tp real(8),dimension(:),allocatable :: l_egy_r,l_egy_tp real(8),dimension(:),allocatable :: m_egy_cmb,l_egy_cmb,m_egy_surface,l_egy_surface real(8),dimension(:,:,:,:),allocatable :: vint,vinp real(8),dimension(:,:,:,:,:),allocatable :: dvint,dvinp real(8),dimension(:),allocatable :: rr,dr,r2inv real(8),dimension(:,:,:),allocatable :: glm,hlm real(8),dimension(:,:),allocatable :: egy1 character(len=8) :: filename3(0:9) character(len=8) :: filename4(10:99) character(len=8) :: filename5(100:127) character(len=9) :: filename30(0:9) character(len=9) :: filename40(10:99) character(len=9) :: filename50(100:127) character(len=8) :: filename300(0:9) character(len=8) :: filename400(10:99) character(len=8) :: filename500(100:127) character(len=9) :: filename3000(0:9) character(len=9) :: filename4000(10:99) character(len=9) :: filename5000(100:127)
character(len=1) :: myrank3 ! 0-9 character(len=2) :: myrank4 ! 10-99 character(len=3) :: myrank5 ! 100-nl
allocate(vint(0:1,0:nl+1,0:nl+1,0:nr),vinp(0:1,0:nl+1,0:nl+1,0:nr) ) allocate(dvint(0:1,0:nl+1,0:nl+1,0:nr,1:2),dvinp(0:1,0:nl+1,0:nl+1,0:nr,1:2) )
allocate(rr(0:nr),r2inv(0:nr),dr(0:nr)) allocate(mp(0:nl+1),lp(0:nl+1),l(0:nl+1,0:nl+1),m(0:nl+1,0:nl+1)) allocate(glm(0:nl+1,0:nl+1,0:nr),hlm(0:nl+1,0:nl+1,0:nr)) allocate(m_egy_r(0:nl+1),m_egy_tp(0:nl+1)) allocate(l_egy_r(0:nl+1),l_egy_tp(0:nl+1)) allocate(m_egy_cmb(0:nl+1),l_egy_cmb(0:nl+1),egy1(0:nl+1,0:nl+1)) allocate(m_egy_surface(0:nl+1),l_egy_surface(0:nl+1)) ri =eta/(1.0d0-eta) !inner core radius ro =1.0d0/(1.0d0-eta) !outer core radius volume=4.0d0*pi/3.0d0*(ro**3-ri**3) !volume of the fluid core
!r方向のgridの設定 do k=0,nr if(k==0) then rr(k)=ri elseif(k==nr) then rr(k)=ro else rr(k)=0.5d0*(ri+ro-cos(pi*k/nr)) endif
if(k>0.and.k<nr) then dr(k)=(ro-ri)*sin(0.5d0*pi/nr)*sin(0.5d0*pi*(2*k-1)/nr) endif r2inv(k)=1.0d0/(rr(k)*rr(k)) end do
!====== mごとのファイルを設定する ====== do mm=0,9 write(myrank3,"(I1)") mm filename3(mm) = var1//"_m.00"//myrank3 filename30(mm) = "d"//var1//"_m.00"//myrank3 filename300(mm) = var2//"_m.00"//myrank3 filename3000(mm) = "d"//var2//"_m.00"//myrank3 open(mm+1000,file=filename3(mm),status="old") open(mm+1200,file=filename30(mm),status="old") open(mm+1400,file=filename300(mm),status="old") open(mm+1600,file=filename3000(mm),status="old") enddo
do mm=10,99 write(myrank4,"(I2)") mm filename4(mm) = var1//"_m.0"//myrank4 filename40(mm) = "d"//var1//"_m.0"//myrank4 filename400(mm) = var2//"_m.0"//myrank4 filename4000(mm) = "d"//var2//"_m.0"//myrank4 open(mm+1000,file=filename4(mm),status="old") open(mm+1200,file=filename40(mm),status="old") open(mm+1400,file=filename400(mm),status="old") open(mm+1600,file=filename4000(mm),status="old") enddo
do mm=100,nl write(myrank5,"(I3)") mm filename5(mm) = var1//"_m."//myrank5 filename50(mm) = "d"//var1//"_m."//myrank5 filename500(mm) = var2//"_m."//myrank5 filename5000(mm) = "d"//var2//"_m."//myrank5 open(mm+1000,file=filename5(mm),status="old") open(mm+1200,file=filename50(mm),status="old") open(mm+1400,file=filename500(mm),status="old") open(mm+1600,file=filename5000(mm),status="old") enddo
!dataを読み込む !========================================================================================= do mm=0,nl do ll=mm,nl do k=0,nr read(mm+1000,*) l(ll,mm),m(ll,mm),kk,vinp(0,ll,mm,k),vinp(1,ll,mm,k) read(mm+1200,*) l(ll,mm),m(ll,mm),kk,dvinp(0,ll,mm,k,1),dvinp(1,ll,mm,k,1),&& dvinp(0,ll,mm,k,2),dvinp(1,ll,mm,k,2)
read(mm+1400,*) l(ll,mm),m(ll,mm),kk,vint(0,ll,mm,k),vint(1,ll,mm,k) read(mm+1600,*) l(ll,mm),m(ll,mm),kk,dvint(0,ll,mm,k,1),dvint(1,ll,mm,k,1),&& dvint(0,ll,mm,k,2),dvint(1,ll,mm,k,2) enddo enddo enddo
!コア平均したエネルギー do mm=0,nl do ll=mm,nl
if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr-1 egy_r = egy_r & & + 0.5d0*dble(l(ll,mm))*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *( ( vinp(0,ll,mm,k)*vinp(0,ll,mm,k)*r2inv(k) + vinp(0,ll,mm,k+1)*vinp(0,ll,mm,k+1)*r2inv(k+1) )&& +( vinp(1,ll,mm,k)*vinp(1,ll,mm,k)*r2inv(k) + vinp(1,ll,mm,k+1)*vinp(1,ll,mm,k+1)*r2inv(k+1) ))
egy_t = egy_t && + 0.5d0*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *((dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(0,ll,mm,k+1,1)*dvinp(0,ll,mm,k+1,1) ) && +(dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1) + dvinp(1,ll,mm,k+1,1)*dvinp(1,ll,mm,k+1,1)) )
egy_p = egy_p && + 0.5d0*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *( ( vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(0,ll,mm,k+1)*vint(0,ll,mm,k+1) ) && +( vint(1,ll,mm,k)*vint(1,ll,mm,k) + vint(1,ll,mm,k+1)*vint(1,ll,mm,k+1) ) ) enddo if(m(ll,mm)==0.and.l(ll,mm)==1) then dip = (egy_r+egy_t) else if(m(ll,mm)==0.and.l(ll,mm)==2) then qua = (egy_r+egy_t)-dip else if(m(ll,mm)==0.and.l(ll,mm)==3) then oct =(egy_r+egy_t)-qua-dip end if enddo enddo
write(*,*) "Energy_Density=", (egy_r + egy_t + egy_p)/(2.0d0*ek*pm*volume) write(*,*) "Dip=", dip/(2.0d0*ek*pm*volume) write(*,*) "Qua=", qua/(2.0d0*ek*pm*volume) write(*,*) "Octa=",oct/(2.0d0*ek*pm*volume)
!エネルギーvs.m を計算する!========================================================================================= open(9,file=var//'egy_m.dat',status='replace')
do mm=0,nl do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr-1 m_egy_r(mm) = m_egy_r(mm) && + (0.5d0*dble(l(ll,mm))*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *( ( vinp(0,ll,mm,k)*vinp(0,ll,mm,k)*r2inv(k) + vinp(0,ll,mm,k+1)*vinp(0,ll,mm,k+1)*r2inv(k+1) )&& +( vinp(1,ll,mm,k)*vinp(1,ll,mm,k)*r2inv(k) + vinp(1,ll,mm,k+1)*vinp(1,ll,mm,k+1)*r2inv(k+1) )))
m_egy_tp(mm) = m_egy_tp(mm) && +(0.5d0*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *((( dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(0,ll,mm,k+1,1)*dvinp(0,ll,mm,k+1,1) ) && +(dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1) + dvinp(1,ll,mm,k+1,1)*dvinp(1,ll,mm,k+1,1)) ) & & + ( ( vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(0,ll,mm,k+1)*vint(0,ll,mm,k+1) ) && +( vint(1,ll,mm,k)*vint(1,ll,mm,k) + vint(1,ll,mm,k+1)*vint(1,ll,mm,k+1) ) ) )) mp(mm)=m(ll,mm) enddo enddo write(9,*) mp(mm),(m_egy_r(mm)+m_egy_tp(mm))/(2.0d0*ek*pm*volume) !m,ポロイダル，トロイダル，全体 enddo
close(9)!=========================================================================================
!エネルギーvs.l を計算する!========================================================================================= open(10,file=var//'egy_l.dat',status='replace')
do mm=0,nl l_egy_r(ll)=0.0d0 l_egy_tp(ll)=0.0d0 do ll=mm,nl
if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif
do k=0,nr-1 l_egy_r(ll) = l_egy_r(ll) && + (0.5d0*dble(l(ll,mm))*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *( ( vinp(0,ll,mm,k)*vinp(0,ll,mm,k)*r2inv(k) + vinp(0,ll,mm,k+1)*vinp(0,ll,mm,k+1)*r2inv(k+1) )&& +( vinp(1,ll,mm,k)*vinp(1,ll,mm,k)*r2inv(k) + vinp(1,ll,mm,k+1)*vinp(1,ll,mm,k+1)*r2inv(k+1) )))
l_egy_tp(ll) = l_egy_tp(ll) && +(0.5d0*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm*dr(k+1) && *((( dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(0,ll,mm,k+1,1)*dvinp(0,ll,mm,k+1,1) ) && +(dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1) + dvinp(1,ll,mm,k+1,1)*dvinp(1,ll,mm,k+1,1)) ) & & + ( ( vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(0,ll,mm,k+1)*vint(0,ll,mm,k+1) ) && +( vint(1,ll,mm,k)*vint(1,ll,mm,k) + vint(1,ll,mm,k+1)*vint(1,ll,mm,k+1) ) ) )) lp(ll)=l(ll,mm) enddo enddo enddo
do ll=0,nl write(10,*) lp(ll),(l_egy_r(ll)+l_egy_tp(ll))/(2.0d0*ek*pm*volume) enddo
close(10)
open(11,file=var//'egy_l_cmb.dat',status='replace') open(12,file=var//'egy_m_cmb.dat',status='replace')
! Lowes spectrum!at CMB do mm=0,nl do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then egy_cmb = egy_cmb && +((dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro*ro*ro) && * (vinp(0,ll,mm,k)*vinp(0,ll,mm,k) + vinp(1,ll,mm,k)*vinp(1,ll,mm,k)))*ro*ro &
& +(dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro) && * ((dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1)) && +(vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(1,ll,mm,k)*vint(1,ll,mm,k))))*ro*ro ) endif enddo enddo enddo write(*,*) "B_average at CMB=", dsqrt(egy_cmb/(2.0d0*ek*pm*4.0d0*pi*ro*ro))
do mm=0,nl do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then m_egy_cmb(mm) = m_egy_cmb(mm) && +((dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro*ro*ro) && * (vinp(0,ll,mm,k)*vinp(0,ll,mm,k) + vinp(1,ll,mm,k)*vinp(1,ll,mm,k)))*ro*ro &
& +(dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro) && * ((dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1)) && +(vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(1,ll,mm,k)*vint(1,ll,mm,k))))*ro*ro ) endif enddo enddo write(12,*) mp(mm),m_egy_cmb(mm)/(2.0d0*ek*pm*4.0d0*pi*ro*ro) enddo
do mm=0,nl l_egy_cmb(ll)=0.0d0 do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then l_egy_cmb(ll) = l_egy_cmb(ll) && +((dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro*ro*ro) && * (vinp(0,ll,mm,k)*vinp(0,ll,mm,k) + vinp(1,ll,mm,k)*vinp(1,ll,mm,k)))*ro*ro &
& +(dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*nm/(ro*ro) && * ((dvinp(0,ll,mm,k,1)*dvinp(0,ll,mm,k,1) + dvinp(1,ll,mm,k,1)*dvinp(1,ll,mm,k,1)) && +(vint(0,ll,mm,k)*vint(0,ll,mm,k) + vint(1,ll,mm,k)*vint(1,ll,mm,k))))*ro*ro ) lp(ll)=l(ll,mm) endif enddo enddo enddo
do ll=0,nl write(11,*) lp(ll),l_egy_cmb(ll)/(2.0d0*ek*pm*4.0d0*pi*ro*ro) enddo
close(11) close(12)
open(13,file=var//'egy_m_surface.dat',status='replace') open(14,file=var//'egy_l_surface.dat',status='replace')
!at surface
!ガウス係数を計算する!========================================= do mm=0,nl do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then glm(ll,mm,k) = dble(l(ll,mm))/(ro*ro) * ((r_core/r_moon)**(dble(l(ll,mm))+2.0d0)) && * dsqrt((2.0d0*dble(l(ll,mm))+1.0d0)/(2.0d0*(2.0d0-delta))) * vinp(0,ll,mm,k)
hlm(ll,mm,k) = dble(l(ll,mm))/(ro*ro) * ((r_core/r_moon)**(dble(l(ll,mm))+2.0d0)) && * dsqrt((2.0d0*dble(l(ll,mm))+1.0d0)/(2.0d0*(2.0d0-delta))) * vinp(1,ll,mm,k) endif enddo enddo enddo
do mm=0,nl do ll=mm,nl
if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then egy_surface = egy_surface && +( (2.0d0*(2.0d0-delta)/(2.0d0*dble(l(ll,mm))+1.0d0))*nm && * (dble(l(ll,mm)) + 1.0d0)*( 2.0d0*dble(l(ll,mm)) + 1.0d0 ) && * (glm(ll,mm,k)*glm(ll,mm,k) + hlm(ll,mm,k)*hlm(ll,mm,k)))*(r_moon*r_moon) endif enddo enddo enddo write(*,*) "B_average at surface=", dsqrt(egy_surface/(2.0d0*ek*pm*4.0d0*pi*r_moon*r_moon))
do mm=0,nl do ll=mm,nl
if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then m_egy_surface(mm) = m_egy_surface(mm) && +( (2.0d0*(2.0d0-delta)/(2.0d0*dble(l(ll,mm))+1.0d0))*nm && * (dble(l(ll,mm)) + 1.0d0)*( 2.0d0*dble(l(ll,mm)) + 1.0d0 ) && * (glm(ll,mm,k)*glm(ll,mm,k) + hlm(ll,mm,k)*hlm(ll,mm,k)))*(r_moon*r_moon) endif enddo enddo write(13,*) mp(mm),m_egy_surface(mm)/(2.0d0*ek*pm*4.0d0*pi*r_moon*r_moon) enddo
do mm=0,nl l_egy_cmb(ll)=0.0d0 do ll=mm,nl if(m(ll,mm)==0) then delta=1 nm=2.0d0*pi else delta=0 nm=pi endif do k=0,nr if(k==nr) then
l_egy_surface(ll) = l_egy_surface(ll) && +( (2.0d0*(2.0d0-delta)/(2.0d0*dble(l(ll,mm))+1.0d0))*nm && * (dble(l(ll,mm)) + 1.0d0)*( 2.0d0*dble(l(ll,mm)) + 1.0d0 ) && * (glm(ll,mm,k)*glm(ll,mm,k) + hlm(ll,mm,k)*hlm(ll,mm,k)))*(r_moon*r_moon) lp(ll)=l(ll,mm) endif enddo enddo enddo
do ll=0,nl write(14,*) lp(ll),l_egy_surface(ll)/(2.0d0*ek*pm*4.0d0*pi*r_moon*r_moon) enddo
close(13) close(14)
do mm=0,nl close(mm+1000) close(mm+1200) close(mm+1400) close(mm+1600) enddo end program egy