Deep Earth Physics - 赤道断面の渦度と磁場のｚ成分を計算するプログラム

赤道断面の渦度と磁場のｚ成分を計算するプログラム

Bz@赤道

program bz_e
implicit none
real(8),parameter :: eta=0.5d0 ! 半径比 !変更部分
character(len=2) :: var1="bp" ! bp,bt,vp,vt,tp,cm !変更部 character(len=2) :: var2="bt" ! bp,bt,vp,vt,tp,cm !変更部 character(len=1) :: varr="b"
integer :: nr=99,nl=127,nt=192,np=400,ph !98:赤道，49:45度断面 integer :: i,j,k,kk,mm,ll,th real(8) :: kaijo,pre_kaijo,sinm,sint,cost,root,time,rd,vmax,vmaxx,vmin,ro,ri real(8), dimension(:,:,:,:),allocatable :: vint real(8), dimension(:,:,:,:,:),allocatable :: dvint real(8), dimension(:,:,:,:),allocatable :: vinp real(8), dimension(:,:,:,:,:),allocatable :: dvinp real(8), dimension(:),allocatable :: r,theta,pt,pp,pr,rinv,invsin,phi real(8), dimension(:,:,:,:),allocatable :: invlt1,invlt2 real(8), dimension(:,:,:),allocatable :: fou1,fou2,plm integer, dimension(:,:),allocatable :: l,m integer,dimension(:),allocatable :: mp1,mp2 real(8), dimension(:,:,:),allocatable :: vout,b_z,b_r,b_theta real(8) :: pi=3.14159265358979323846264338327950288d0
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) :: myrank1 ! 0-9 character(len=2) :: myrank2 ! 10-31 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) ) allocate(dvint(0:1,0:nl+1,0:nl+1,0:nr,1:2) ) allocate(vinp(0:1,0:nl+1,0:nl+1,0:nr) ) allocate(dvinp(0:1,0:nl+1,0:nl+1,0:nr,1:2) )
allocate(r(0:nr),rinv(0:nr),theta(0:nt-1),phi(0:np-1)) allocate(pr(0:nr),pt(0:nt-1),pp(0:np-1),invsin(0:nt-1)) allocate(fou1(0:nr,0:nt-1,0:np-1),fou2(0:nr,0:nt-1,0:np-1),vout(0:nr,0:nt-1,0:np-1)) allocate(b_z(0:nr,0:nt-1,0:np-1),b_r(0:nr,0:nt-1,0:np-1),b_theta(0:nr,0:nt-1,0:np-1)) allocate(invlt1(0:nr,0:nt-1,0:1,0:nl+1),invlt2(0:nr,0:nt-1,0:1,0:nl+1)) allocate(plm(0:nt-1,0:nl+1,0:nl+1)) allocate(mp2(0:nl+1),mp1(0:nl+1),l(0:nl+1,0:nl+1),m(0:nl+1,0:nl+1))
!ラジアンから度に変更する（描画用） rd=180.0d0/pi
!おえかきするときのrのグリッド do k=0,nr ri=eta/(1.0d0-eta) ro=1.0d0/(1.0d0-eta) if(k==0) then r(k)=ri elseif(k==nr) then r(k)=ro else r(k)=0.5d0*(ri+ro-dcos(pi*dble(k)/dble(nr))) endif rinv(k)=1.0d0/r(k) enddo
!(sinθ)^-1 の計算 do th=0,nt-1 if(th==0.and.th==nt-1) then invsin(th)=0.0d0 else invsin(th)=1.0d0/dsin(pi*dble(th)/dble(nt-1)) endif enddo!!!====== mごとのファイルを設定する ====== if(nl<100) then 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,nl 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
elseif(nl>=100) then 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 endif

!===== 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 if(m(ll,mm)>nl) then write(*,*) 'm is larger than nl' !mが打ち切り次数nlより大きいとダメ!! stop end if do ll=mm,nl
kaijo=dsqrt(0.5d0) do j=1,mm kaijo=kaijo*dsqrt(dble(2*j+1)/dble(2*j)) enddo
!l=mのとき---------------------------------------------------! if( l(ll,mm) == m(ll,mm) ) then do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) sinm = sint**m(ll,mm) plm(th,ll,mm) = kaijo * sinm enddo !------------------------------------------------------------!! !l=m+1のとき-------------------------------------------------! elseif( l(ll,mm)==(m(ll,mm)+1) ) then do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) root = dsqrt((2.0d0*dble(m(ll,mm)) + 3.0d0)) plm(th,ll,mm) = cost*root*plm(th,ll-1,mm) enddo !------------------------------------------------------------!!! !それ以外(ルジャンドル陪関数の一般解)---------------------------------------------------! else do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) plm(th,ll,mm) = sqrt((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0) && / ((dble(l(ll,mm))+dble(m(ll,mm)))*(dble(l(ll,mm))-dble(m(ll,mm))))) && *cost*plm(th,ll-1,mm) &
& - sqrt((dble(l(ll,mm))+dble(m(ll,mm))-1.0d0) && *(dble(l(ll,mm))-dble(m(ll,mm))-1.0d0)*(2.0d0*dble(l(ll,mm))+1.0d0) && / ((dble(l(ll,mm))+dble(m(ll,mm))) && *(dble(l(ll,mm))-dble(m(ll,mm)))*(2.0d0*dble(l(ll,mm))-3.0d0))) && *plm(th,ll-2,mm) enddo endif !---------------------------------------------------------------------------------------! enddo enddo
!-----------------------------------------------------------------------------------------------!!-------------------------------------b_r成分の計算---------------------------------------------!!-----------------------------------------------------------------------------------------------!!===== mごとにルジャンドル逆変換を実行する ================ do th=0,nt-1 do k=0,nr do mm=0,nl do ll=mm,nl invlt1(k,th,0,mm) = invlt1(k,th,0,mm) && + (dble(l(ll,mm)*(dble(l(ll,mm))+1.0d0))*vinp(0,ll,mm,k)*plm(th,ll,mm)) invlt1(k,th,1,mm) = invlt1(k,th,1,mm) && + (dble(l(ll,mm)*(dble(l(ll,mm))+1.0d0))*vinp(1,ll,mm,k)*plm(th,ll,mm))
mp1(mm)=m(ll,mm) enddo enddo enddo enddo!!!===== フーリエ変換の和を計算する ================ do ph=0,np-1 phi=2.0d0*pi*dble(ph)/dble(np-1) do th=0,nt-1 do k=0,nr do mm=0,nl fou1(k,th,ph) = fou1(k,th,ph) && +( invlt1(k,th,0,mm)*dcos(mp1(mm)*phi(ph)) + invlt1(k,th,1,mm)*dsin(mp1(mm)*phi(ph)) ) enddo enddo enddo enddo!!===== 物理空間の変数に変換する ========================! do ph=0,np-1 do th=0,nt-1 do k=0,nr ! b_r(k,th,ph)=(rinv(k)*rinv(k)*fou1(k,th,ph))*dcos(-dble(th)*pi/dble(nt-1)+0.5d0*pi) b_r(k,th,ph)=rinv(k)*rinv(k)*fou1(k,th,ph) enddo enddo enddo

!-----------------------------------------------------------------------------------------------!!-------------------------------------b_theta成分の計算-----------------------------------------!!-----------------------------------------------------------------------------------------------!!===== mごとにルジャンドル逆変換を実行する ================ do th=0,nt-1 do k=0,nr do mm=0,nl do ll=mm,nl invlt2(k,th,0,mm) = invlt2(k,th,0,mm) + &
& ( (((l(ll,mm)-1.0d0)*dsqrt((dble(l(ll,mm))+dble(m(ll,mm)))*(dble(l(ll,mm))-dble(m(ll,mm))) && /((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0)))*(dvinp(0,ll-1,mm,k,1))) &
& -( (l(ll,mm)+2.0d0)*dsqrt((dble(l(ll+1,mm))+dble(m(ll+1,mm)))*(dble(l(ll+1,mm))-dble(m(ll+1,mm))) && /((2.0d0*dble(l(ll+1,mm))+1.0d0)*(2.0d0*dble(l(ll+1,mm))-1.0d0)))*(dvinp(0,ll+1,mm,k,1)) ) &
& + (m(ll,mm)*vint(1,ll,mm,k))) * plm(th,ll,mm) )

invlt2(k,th,1,mm) = invlt2(k,th,1,mm) + && ( (((l(ll,mm)-1.0d0)*dsqrt((dble(l(ll,mm))+1.0d0)*(dble(l(ll,mm))-dble(m(ll,mm))) && /((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0)))*(dvinp(1,ll-1,mm,k,1))) &
& -( (l(ll,mm)+2.0d0)*dsqrt((dble(l(ll+1,mm))+dble(m(ll+1,mm)))*(dble(l(ll+1,mm))-dble(m(ll+1,mm))) && /((2.0d0*dble(l(ll+1,mm))+1.0d0)*(2.0d0*dble(l(ll+1,mm))-1.0d0)))*(dvinp(1,ll+1,mm,k,1)) ) &
& - (m(ll,mm)*vint(0,ll,mm,k))) * plm(th,ll,mm) ) mp2(mm)=m(ll,mm) enddo enddo enddo enddo

!===== フーリエ変換の和を計算する ================ do ph=0,np-1 phi=2.0d0*pi*dble(ph)/dble(np-1) do th=0,nt-1 do k=0,nr do mm=0,nl fou2(k,th,ph) = fou2(k,th,ph) && +( invlt2(k,th,0,mm)*dcos(mp2(mm)*phi(ph)) + invlt2(k,th,1,mm)*dsin(mp2(mm)*phi(ph)) ) enddo enddo enddo enddo
!===== 物理空間の変数に変換する ========================! do ph=0,np-1 do th=0,nt-1 do k=0,nr b_theta(k,th,ph)=rinv(k)*fou2(k,th,ph)*invsin(th) enddo enddo enddo
!-----------------------------------------------------------------------------------------------!!-------------------------------------b_z成分の計算---------------------------------------------!!-----------------------------------------------------------------------------------------------! do ph=0,np-1 do th=0,nt-1 do k=0,nr theta(th)=pi*dble(th)/dble(nt-1) b_z(k,th,ph)=b_r(k,th,ph)*dcos(theta(th)) - b_theta(k,th,ph)*dsin(theta(th)) enddo enddo enddo
do ph=0,np-1 do th=0,nt-1 do k=0,nr if(th==96) vout(k,th,ph)=(b_z(k,th,ph)) enddo enddo enddo
!最大値と最小値を検索 do ph=0,np-1 do th=0,nt-1 if(th==96) then do k=0,nr if( vmax < vout(k,th,ph) ) then vmax = vout(k,th,ph) elseif( vmax > vout(k,th,ph) ) then vmax=vmax endif enddo endif enddo enddo write(*,*) "vmax=",vmax
do ph=0,np-1 do th=0,nt-1 if(th==96) then do k=0,nr if( vmin > vout(k,th,ph) ) then vmin = vout(k,th,ph) elseif( vmin < vout(k,th,ph) ) then vmin=vmin endif enddo endif enddo enddo write(*,*) "vmin=",vmin

!結果の出力 ++++++++++++++++++++++++++++++++++++++++++++ open(8,file=varr//"_z.dat",status="replace")
do ph=0,np-1 if( dble(ph)<=dble((np-1)/2) ) then pp(ph)=2.0d0*pi*dble(ph)/dble(np-1) else pp(ph)=-2.0d0*pi*dble((np-1)-ph)/dble(np-1) endif
do th=0,nt-1 if(th==96) then do k=0,nr pr(k)=(1.0d1 * ( 1.0d0 - eta ) * r(k)) if(vmax>dabs(vmin)) then write(8,"(3d13.5)") pp(ph)*rd,pr(k),vout(k,th,ph)/vmax else write(8,"(3d13.5)") pp(ph)*rd,pr(k),vout(k,th,ph)/dabs(vmin) endif enddo endif enddo enddo
close(8)

open(2,file="maxmin.dat",status="unknown",position="append") write(2,*) varr,"z_component at equatorial plane:","max=",vmax,"min=",vmin close(2)
end program bz_e

omgz@赤道

program conversion_r
implicit none
real(8),parameter :: eta=0.5d0 ! 半径比 !変更部分 character(len=2) :: var1="vp" ! bp,bt,vp,vt,tp,cm !変更部 character(len=2) :: var2="vt" ! bp,bt,vp,vt,tp,cm !変更部 character(len=1) :: var="v" integer :: nr=99,nl=127,nt=192,np=400,ph integer :: i,j,k,kk,mm,ll,th real(8) :: kaijo,sinm,sint,cost,root,time,rd,vmax,vmaxx,vmin,ro,ri
real(8), dimension(:,:,:,:),allocatable :: vint real(8), dimension(:,:,:,:,:),allocatable :: dvint real(8), dimension(:,:,:,:),allocatable :: vinp real(8), dimension(:,:,:,:,:),allocatable :: dvinp real(8), dimension(:),allocatable :: r,theta,pt,pp,pr,rinv,invsin,phi real(8), dimension(:,:,:,:),allocatable :: invlt1,invlt2 real(8), dimension(:,:,:),allocatable :: fou1,fou2,plm integer, dimension(:,:),allocatable :: l,m integer,dimension(:),allocatable :: mp1,mp2 real(8), dimension(:,:,:),allocatable :: omgz,omgr,omgt,vout real(8) :: pi=3.14159265358979323846264338327950288d0
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) ) allocate(dvint(0:1,0:nl+1,0:nl+1,0:nr,1:2) ) allocate(vinp(0:1,0:nl+1,0:nl+1,0:nr) ) allocate(dvinp(0:1,0:nl+1,0:nl+1,0:nr,1:2) )
allocate(r(0:nr),rinv(0:nr),theta(0:nt-1),phi(0:np-1)) allocate(pr(0:nr),pt(0:nt-1),pp(0:np-1),invsin(0:nt-1)) allocate(fou1(0:nr,0:nt-1,0:np-1),fou2(0:nr,0:nt-1,0:np-1)) allocate(omgt(0:nr,0:nt-1,0:np-1)) allocate(omgr(0:nr,0:nt-1,0:np-1)) allocate(omgz(0:nr,0:nt-1,0:np-1)) allocate(vout(0:nr,0:nt-1,0:np-1)) allocate(invlt1(0:nr,0:nt-1,0:1,0:nl+1),invlt2(0:nr,0:nt-1,0:1,0:nl+1)) allocate(plm(0:nt-1,0:nl+1,0:nl+1)) allocate(mp2(0:nl+1),mp1(0:nl+1),l(0:nl+1,0:nl+1),m(0:nl+1,0:nl+1))
!ラジアンから度に変更する（描画用） rd=180.0d0/pi
!おえかきするときのrのグリッド do k=0,nr ri=eta/(1.0d0-eta) ro=1.0d0/(1.0d0-eta) if(k==0) then r(k)=ri elseif(k==nr) then r(k)=ro else r(k)=0.5d0*(ri+ro-dcos(pi*dble(k)/dble(nr))) endif rinv(k)=1.0d0/r(k) enddo
!(sinθ)^-1 の計算 do th=0,nt-1 if(th==0.and.th==nt-1) then invsin(th)=0.0d0 else invsin(th)=1.0d0/dsin(pi*dble(th)/dble(nt-1)) endif enddo!!!====== mごとのファイルを設定する ====== if(nl<100) then 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,nl 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
elseif(nl>=100) then 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 endif

!===== 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 if(m(ll,mm)>nl) then write(*,*) 'm is larger than nl' !mが打ち切り次数nlより大きいとダメ!! stop end if do ll=mm,nl
kaijo=dsqrt(0.5d0) do j=1,mm kaijo=kaijo*dsqrt(dble(2*j+1)/dble(2*j)) enddo
!l=mのとき---------------------------------------------------! if( l(ll,mm) == m(ll,mm) ) then do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) sinm = sint**m(ll,mm) plm(th,ll,mm) = kaijo * sinm enddo !------------------------------------------------------------!! !l=m+1のとき-------------------------------------------------! elseif( l(ll,mm)==(m(ll,mm)+1) ) then do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) root = dsqrt((2.0d0*dble(m(ll,mm)) + 3.0d0)) plm(th,ll,mm) = cost*root*plm(th,ll-1,mm) enddo !------------------------------------------------------------!!! !それ以外(ルジャンドル陪関数の一般解)---------------------------------------------------! else do th=0,nt-1 theta(th)=pi*dble(th)/dble(nt-1) cost = dcos(theta(th)) sint = dsin(theta(th)) plm(th,ll,mm) = sqrt((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0) && / ((dble(l(ll,mm))+dble(m(ll,mm)))*(dble(l(ll,mm))-dble(m(ll,mm))))) && *cost*plm(th,ll-1,mm) &
& - sqrt((dble(l(ll,mm))+dble(m(ll,mm))-1.0d0) && *(dble(l(ll,mm))-dble(m(ll,mm))-1.0d0)*(2.0d0*dble(l(ll,mm))+1.0d0) && / ((dble(l(ll,mm))+dble(m(ll,mm))) && *(dble(l(ll,mm))-dble(m(ll,mm)))*(2.0d0*dble(l(ll,mm))-3.0d0))) && *plm(th,ll-2,mm) enddo endif !---------------------------------------------------------------------------------------! enddo enddo
!-----------------------------------------------------------------------------------------------!!-------------------------------------omg_r成分の計算-------------------------------------------!!-----------------------------------------------------------------------------------------------!!===== mごとにルジャンドル逆変換を実行する ================ do th=0,nt-1 do k=0,nr do mm=0,nl do ll=mm,nl invlt1(k,th,0,mm) = invlt1(k,th,0,mm) && + (dble(l(ll,mm)*(dble(l(ll,mm))+1.0d0))*vint(0,ll,mm,k)*plm(th,ll,mm)) invlt1(k,th,1,mm) = invlt1(k,th,1,mm) && + (dble(l(ll,mm)*(dble(l(ll,mm))+1.0d0))*vint(1,ll,mm,k)*plm(th,ll,mm)) mp1(mm)=m(ll,mm) enddo enddo enddo enddo!!===== フーリエ変換の和を計算する =============== do ph=0,np-1 phi=2.0d0*pi*dble(ph)/dble(np-1) do th=0,nt-1 do k=0,nr do mm=0,nl fou1(k,th,ph) = fou1(k,th,ph) && +( invlt1(k,th,0,mm)*dcos(mp1(mm)*phi(ph)) + invlt1(k,th,1,mm)*dsin(mp1(mm)*phi(ph)) ) enddo enddo enddo enddo
!===== 物理空間の変数に変換する ========================! do ph=0,np-1 do th=0,nt-1 do k=0,nr omgr(k,th,ph)=rinv(k)*rinv(k)*fou1(k,th,ph) enddo enddo enddo! !!!-----------------------------------------------------------------------------------------------!!-------------------------------------omg_theta成分の計算---------------------------------------!!-----------------------------------------------------------------------------------------------!!===== mごとにルジャンドル逆変換を実行する ================ do th=0,nt-1 do k=0,nr do mm=0,nl do ll=mm,nl invlt2(k,th,0,mm) = invlt2(k,th,0,mm) + &
& ( (((l(ll,mm)-1.0d0)*dsqrt((dble(l(ll,mm))+dble(m(ll,mm)))*(dble(l(ll,mm))-dble(m(ll,mm))) && /((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0)))*(dvint(0,ll-1,mm,k,1))) &
& -( (l(ll,mm)+2.0d0)*dsqrt((dble(l(ll+1,mm))+dble(m(ll+1,mm)))*(dble(l(ll+1,mm))-dble(m(ll+1,mm)))&& /((2.0d0*dble(l(ll+1,mm))+1.0d0)*(2.0d0*dble(l(ll+1,mm))-1.0d0)))*(dvint(0,ll+1,mm,k,1)) ) &
& -(m(ll,mm) && *(dvinp(1,ll+1,mm,k,2)-(dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*rinv(k)*rinv(k)*vinp(1,ll,mm,k)))))& * plm(th,ll,mm) )

invlt2(k,th,1,mm) = invlt2(k,th,1,mm) + & & ( (((l(ll,mm)-1.0d0)*dsqrt((dble(l(ll,mm))+dble(m(ll,mm)))*(dble(l(ll,mm))-dble(m(ll,mm))) && /((2.0d0*dble(l(ll,mm))+1.0d0)*(2.0d0*dble(l(ll,mm))-1.0d0)))*(dvinp(1,ll-1,mm,k,1))) &
& -( (l(ll,mm)+2.0d0)*dsqrt((dble(l(ll+1,mm))+dble(m(ll+1,mm)))*(dble(l(ll+1,mm))-dble(m(ll+1,mm)))&& /((2.0d0*dble(l(ll+1,mm))+1.0d0)*(2.0d0*dble(l(ll+1,mm))-1.0d0)))*(dvinp(1,ll+1,mm,k,1)) ) &
& +(m(ll,mm) && *(dvinp(0,ll+1,mm,k,2)-(dble(l(ll,mm))*(dble(l(ll,mm))+1.0d0)*rinv(k)*rinv(k)*vinp(0,ll,mm,k)))))& * plm(th,ll,mm) ) mp2(mm)=m(ll,mm) enddo enddo enddo enddo
!===== フーリエ変換の和を計算する ================ do ph=0,np-1 phi=2.0d0*pi*dble(ph)/dble(np-1) do th=0,nt-1 do k=0,nr do mm=0,nl fou2(k,th,ph) = fou2(k,th,ph) && +( invlt2(k,th,0,mm)*dcos(mp2(mm)*phi(ph)) + invlt2(k,th,1,mm)*dsin(mp2(mm)*phi(ph)) ) enddo enddo enddo enddo
!===== 物理空間の変数に変換する ========================! do ph=0,np-1 do th=0,nt-1 do k=0,nr omgt(k,th,ph)=rinv(k)*fou2(k,th,ph)*invsin(th) enddo enddo enddo

!-----------------------------------------------------------------------------------------------!!-------------------------------------omg_z成分の計算-------------------------------------------!!-----------------------------------------------------------------------------------------------! do ph=0,np-1 do th=0,nt-1 do k=0,nr theta(th)=pi*dble(th)/dble(nt-1) omgz(k,th,ph)=omgr(k,th,ph)*dcos(theta(th))-omgt(k,th,ph)*dsin(theta(th)) enddo enddo enddo
do ph=0,np-1 do th=0,nt-1 do k=0,nr if(th==96) vout(k,th,ph)=omgz(k,th,ph) enddo enddo enddo
!最大値と最小値を検索 do ph=0,np-1 do th=0,nt-1 if(th==96) then do k=0,nr if( vmax < vout(k,th,ph) ) then vmax = vout(k,th,ph) elseif( vmax > vout(k,th,ph) ) then vmax=vmax endif enddo endif enddo enddo write(*,*) "vmax=",vmax
do ph=0,np-1 do th=0,nt-1 if(th==96) then do k=0,nr if( vmin > vout(k,th,ph) ) then vmin = vout(k,th,ph) elseif( vmin < vout(k,th,ph) ) then vmin=vmin endif enddo endif enddo enddo write(*,*) "vmin=",vmin

!結果の出力 ++++++++++++++++++++++++++++++++++++++++++++
open(8,file="omgz_e.dat",status="replace") do ph=0,np-1 if( dble(ph)<=dble((np-1)/2) ) then pp(ph)=2.0d0*pi*dble(ph)/dble(np-1) else pp(ph)=-2.0d0*pi*dble((np-1)-ph)/dble(np-1) endif
do th=0,nt-1 if(th==96) then do k=0,nr pr(k)=(1.0d1 * ( 1.0d0 - eta ) * r(k)) if(vmax>dabs(vmin)) then write(8,"(3d13.5)") pp(ph)*rd,pr(k), vout(k,th,ph)/vmax else write(8,"(3d13.5)") pp(ph)*rd,pr(k), vout(k,th,ph)/dabs(vmin) endif enddo endif enddo enddo close(8)

open(2,file="maxmin.dat",status="unknown",position="append") write(2,*) "omgz_component at meriodial plane:","max=",vmax,"min=",vmin close(2)
end program conversion_r

ＧＭＴでお絵かき

(Ver:GMT5.4.1です．GMT6では動かない（たぶん）)

#! /bin/csh
rm z.ps
set input1 = omgz_e.datset input2 = b_z.dat
set output = z.ps
set range = -180./180./5./10.
set tmp1 = omg_z.grdset tmp2 = b_z.grd
set fcolor = clr.cpt
gmtset MEASURE_UNIT cm DOTS_PR_INCH 360
makecpt -Cpolar -T-1.0/1.0/0.03 > $fcolor surface $input1 -G$tmp1 -R$range -C -I0.05/0.01 -T0.35surface $input2 -G$tmp2 -R$range -C -I0.05/0.01 -T0.35
grdimage $tmp1 -JP10 -R$range -C$fcolor -X3 -Y3 -P -K >! $output
psscale -C$fcolor -D5c/-1c/10c/1ch -B50 -N320 -O -P -K >> $output
grdcontour $tmp2 -JP10 -R$range -C0.1 -L0.000005/1.0 -W0.9,red -O -P -K >> $outputgrdcontour $tmp2 -JP10 -R$range -C0.1 -L-1.0/-0.000005 -W0.9,blue,- -O -P -K >> $output
psbasemap -JP10 -R$range -Bg60/g2 -O -P >> $output
gmt psconvert -A -Tg z.ps
\rm .gmt* *.conf *.grd #*.cpt

↓こんな図が描ける。（例は半径比0.5のケース）

左図は渦度を色で，，磁場をコンターで書いている。磁力線と重なって渦度が見づらいので，渦度だけの図も描いたりしている（左）。上のお絵かきシェルスクリプトで磁場のお絵かき部分を#でコメントアウトすればよい．

Kageyama and Sato(1997)で言及されているような，負の渦度の所に磁力線が集まっているのが見える。これはRa=6Racrit程度でそんなにレイリー数が大きくないので形が結構きれいにみえている．

色は赤道断面の最大値で規格化してあるので，上のプログラムにあるように，maxmin.datというファイルに出力された最大値を後から下図に付け足さないといけない（powerpointとかで）．あまり効率がよくないが，いちいちシェルスクリプトを書き直すのも面倒なので．．．

これにgrdvectorを使って矢印を付け加えたプログラムも作っているので，そのうち公開する予定．