38 double complex,
allocatable :: flm(:,:),Alm(:,:),Blm(:,:)
39 double precision,
allocatable :: Rlm(:,:), xrg(:)
40 double precision,
public ::
r_s=2.5d0,
r_0=1.d0
41 integer,
allocatable :: lmaxarray(:)
43 logical,
public ::
trunc=.false.
55 double precision,
allocatable :: b_r0(:,:)
56 double precision,
allocatable :: theta(:),phi(:),cfwm(:)
57 double precision :: rsl,xrl,dxr
58 integer :: xm,ym,
l,m,amode,file_handle,il,ir,nlarr,nsh
59 integer,
dimension(MPI_STATUS_SIZE) :: statuss
60 character(len=*) :: mapname
61 character(len=80) :: fharmcoef
64 fharmcoef=mapname//
'.coef'
65 inquire(file=fharmcoef, exist=aexist)
68 call mpi_file_open(mpi_comm_self,fharmcoef,mpi_mode_rdonly, &
69 mpi_info_null,file_handle,
ierrmpi)
70 call mpi_file_read(file_handle,
lmax,1,mpi_integer,statuss,
ierrmpi)
72 call mpi_file_read(file_handle,flm,(
lmax+1)*(
lmax+1),&
73 mpi_double_complex,statuss,
ierrmpi)
74 call mpi_file_close(file_handle,
ierrmpi)
84 inquire(file=mapname,exist=aexist)
86 if(
mype==0)
write(*,
'(2a)')
"can not find file:",mapname
87 call mpistop(
"no input magnetogram found")
89 call mpi_file_open(mpi_comm_self,mapname,mpi_mode_rdonly,mpi_info_null,&
91 call mpi_file_read(file_handle,xm,1,mpi_integer,statuss,
ierrmpi)
92 call mpi_file_read(file_handle,ym,1,mpi_integer,statuss,
ierrmpi)
98 call mpi_file_read(file_handle,theta,ym,mpi_double_precision,&
100 call mpi_file_read(file_handle,phi,xm,mpi_double_precision,&
102 call mpi_file_read(file_handle,b_r0,xm*ym,mpi_double_precision,&
104 call mpi_file_close(file_handle,
ierrmpi)
105 print*,
'nphi,ntheta',xm,ym
106 print*,
'theta range:',minval(theta),maxval(theta)
107 print*,
'phi range:',minval(phi),maxval(phi)
108 print*,
'Brmax,Brmin',maxval(b_r0),minval(b_r0)
112 call coef(b_r0,xm,ym,dcos(theta),dsin(theta),cfwm)
116 amode=ior(mpi_mode_create,mpi_mode_wronly)
117 call mpi_file_open(mpi_comm_self,fharmcoef,amode, &
118 mpi_info_null,file_handle,
ierrmpi)
119 call mpi_file_write(file_handle,
lmax,1,mpi_integer,statuss,
ierrmpi)
120 call mpi_file_write(file_handle,flm,(
lmax+1)*(
lmax+1),&
121 mpi_double_complex,statuss,
ierrmpi)
122 call mpi_file_close(file_handle,
ierrmpi)
128 if(
npe>1)
call mpi_bcast(flm,(
lmax+1)*(
lmax+1),mpi_double_complex,0,&
133 allocate(lmaxarray(nlarr))
137 dxr=(
r_s-
r_0)/dble(nlarr-1)
139 xrg(ir)=dxr*dble(ir-1)+
r_0
158 rlm(
l,m)=dsqrt(dble(
l**2-m**2)/dble(4*
l**2-1))
159 blm(
l,m)=-flm(
l,m)/(1.d0+dble(
l)+dble(
l)*rsl)
160 alm(
l,m)=-rsl*blm(
l,m)
169 integer,
intent(in) :: ym
170 double precision,
intent(in) :: miu(ym)
171 double precision,
intent(out) :: cfwm(ym)
173 double precision,
dimension(ym) :: Pl,Pm2,Pm1,Pprime,sintheta
174 double precision :: lr
177 sintheta=dsqrt(1.d0-miu**2)
184 pl=(2.d0-lr)*pm1*miu-(1.d0-lr)*pm2
189 pprime=(dble(ym)*pl)/sintheta**2
190 cfwm=2.d0/(sintheta*pprime)**2
195 subroutine coef(b_r0,xm,ym,miu,mius,cfwm)
198 integer,
intent(in) :: xm,ym
199 double precision,
intent(in) :: b_r0(xm,ym),cfwm(ym),miu(ym),mius(ym)
201 double complex :: Bm(0:xm-1,0:ym-1)
202 double precision,
dimension(xm) :: fftmr,fftmi
203 double precision,
dimension(0:lmax) :: N_mm
204 double precision,
dimension(ym) :: P_lm1,P_lm2,old_Pmm,P_l
205 double precision :: mr,lr,c1,c2
206 integer :: l,m,i,j,stat
210 fftmr=b_r0(:,i)/dble(xm)
212 call fft(fftmr,fftmi,xm,xm,xm,-1)
213 bm(:,i-1)=(fftmr+(0.d0,1.d0)*fftmi)
215 n_mm(0)=1.d0/dsqrt(4.d0*dpi)
217 n_mm(m)=-n_mm(m-1)*dsqrt(1.d0+1.d0/dble(2*m))
221 p_lm1=p_lm2*miu*dsqrt(3.d0)
223 flm(0,0)=sum(bm(0,:)*p_lm2*cfwm)
225 flm(1,0)=sum(bm(0,:)*p_lm1*cfwm)
228 c1=dsqrt(4.d0-1.d0/lr**2)
229 c2=-(1.d0-1.d0/lr)*dsqrt((2.d0*lr+1.d0)/(2.d0*lr-3.d0))
230 p_l=c1*miu*p_lm1+c2*p_lm2
232 flm(l,0)=sum(bm(0,:)*p_l*cfwm)
243 p_lm2=old_pmm*mius*n_mm(m)/n_mm(m-1)
244 p_lm1=p_lm2*miu*dsqrt(dble(2*m+3))
248 flm(m,m)=sum(bm(m,:)*p_lm2*cfwm)
250 if(m<
lmax) flm(m+1,m)=sum(bm(m,:)*p_lm1*cfwm)
254 c1=dsqrt((4.d0*lr**2-1.d0)/(lr**2-mr**2))
255 c2=-dsqrt(((2.d0*lr+1.d0)*((lr-1.d0)**2-mr**2))/((2.d0*lr-3.d0)*(lr**2-&
257 p_l=c1*miu*p_lm1+c2*p_lm2
258 flm(l,m)=sum(bm(m,:)*p_l*cfwm)
266 subroutine pfss(ixI^L,ixO^L,Bpf,x)
269 integer,
intent(in) :: ixi^
l,ixo^
l
270 double precision,
intent(in) :: x(ixi^s,1:
ndim)
271 double precision,
intent(out) :: bpf(ixi^s,1:
ndir)
273 double complex :: bt(0:
lmax,0:
lmax,ixomin1:ixomax1)
274 double precision :: phase(ixi^s,1:
ndir),bpfiv(ixomin3:ixomax3,ixomin2:ixomax2)
275 double precision :: miu(ixomin2:ixomax2),mius(ixomin2:ixomax2),xr
276 integer ::
l,m,ix^
d,j,l1,l2,ntheta,nphi,ir,qlmax
279 nphi=ixomax3-ixomin3+1
280 ntheta=ixomax2-ixomin2+1
281 miu(ixomin2:ixomax2)=dcos(x(ixomin1,ixomax2:ixomin2:-1,ixomin3,2))
282 mius(ixomin2:ixomax2)=dsin(x(ixomin1,ixomax2:ixomin2:-1,ixomin3,2))
283 do ix1=ixomin1,ixomax1
284 xr=x(ix1,ixomin2,ixomin3,1)
286 do ir=1,
size(lmaxarray)
289 if(ir>
size(lmaxarray)) ir=
size(lmaxarray)
297 bt(
l,m,ix1)=alm(
l,m)*dble(
l)*xr**(
l-1)-blm(
l,m)*dble(
l+1)*xr**(-
l-2)
301 ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
302 do ix3=ixomin3,ixomax3
303 do ix2=ixomin2,ixomax2
304 bpf(ix1,ix2,ix3,1)=bpfiv(ix3,ixomax2-ix2+ixomin2)
311 bt(
l,m,ix1)=-rlm(
l+1,m)*dble(
l+2)*&
312 (alm(
l+1,m)*xr**
l+blm(
l+1,m)*xr**(-
l-3))
313 else if (
l>=1 .and.
l<=
lmax-1)
then
314 bt(
l,m,ix1)=rlm(
l,m)*&
315 dble(
l-1)*(alm(
l-1,m)*xr**(
l-2)+blm(
l-1,m)*&
316 xr**(-
l-1))-rlm(
l+1,m)*dble(
l+2)*&
317 (alm(
l+1,m)*xr**
l+blm(
l+1,m)*xr**(-
l-3))
319 bt(
l,m,ix1)=rlm(
l,m)*&
320 dble(
l-1)*(alm(
l-1,m)*xr**(
l-2)+blm(
l-1,m)*xr**(-
l-1))
325 ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
326 do ix3=ixomin3,ixomax3
327 do ix2=ixomin2,ixomax2
328 bpf(ix1,ix2,ix3,2)=bpfiv(ix3,ixomax2-ix2+ixomin2)/mius(&
336 bt(
l,m,ix1)=(0.d0,1.d0)*m*(alm(
l,m)*xr**(
l-1)+blm(
l,m)*xr**(-
l-2))
340 ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
341 do ix3=ixomin3,ixomax3
342 do ix2=ixomin2,ixomax2
343 bpf(ix1,ix2,ix3,3)=bpfiv(ix3,ixomax2-ix2+ixomin2)/mius(&
363 integer,
intent(in) :: nphi,ntheta,qlmax
364 double complex,
intent(in) :: Bt(0:lmax,0:lmax)
365 double precision,
intent(in) :: phi(nphi),miu(ntheta),mius(ntheta)
366 double precision,
intent(out) :: Bpf(nphi,ntheta)
368 double precision,
dimension(0:lmax,0:lmax) :: cp,phase,Bamp
369 double precision,
dimension(ntheta) :: cp_1_0,cp_l_0,cp_lm1_0,cp_lm2_0,cp_m_m
370 double precision,
dimension(ntheta) :: cp_1_m,cp_l_m,cp_lm1_m,cp_lm2_m,cp_mp1_m
371 double precision :: angpart(nphi)
372 double precision :: ld,md,c1,c2,cp_0_0
373 integer :: l,m,iph,ith
377 phase=atan2(dimag(bt),dble(bt))
381 cp_0_0=dsqrt(1.d0/(4.d0*dpi))
383 bpf=bpf+bamp(0,0)*dcos(phase(0,0))*cp_0_0
386 cp_1_0=dsqrt(3.d0)*miu*cp_0_0
388 bpf(iph,:)=bpf(iph,:)+bamp(1,0)*dcos(phase(1,0))*cp_1_0
398 c1=dsqrt(4.d0*ld**2-1.d0)/ld
399 c2=dsqrt((2.d0*ld+1.d0)/(2.d0*ld-3.d0))*(ld-1.d0)/ld
400 cp_l_0=c1*miu*cp_lm1_0-c2*cp_lm2_0
402 bpf(iph,:)=bpf(iph,:)+bamp(l,0)*dcos(phase(l,0))*cp_l_0
411 cp_m_m=-dsqrt(1.d0+1.d0/(2.d0*md))*mius*cp_m_m
413 angpart(iph)=dcos(md*phi(iph)+phase(m,m))
417 bpf(iph,ith)=bpf(iph,ith)+bamp(m,m)*angpart(iph)*cp_m_m(ith)
423 cp_mp1_m=dsqrt(2.d0*md+3.d0)*miu*cp_m_m
424 angpart=dcos(md*phi+phase(m+1,m))
427 bpf(iph,ith)=bpf(iph,ith)+bamp(m+1,m)*angpart(iph)*cp_mp1_m(ith)
440 c1=dsqrt((4.d0*ld**2-1.d0)/(ld**2-md**2))
441 c2=dsqrt((2.d0*ld+1.d0)*((ld-1.d0)**2-md**2)/(2.d0*ld-3.d0)/(ld**2-md**2))
442 cp_l_m=c1*miu*cp_lm1_m-c2*cp_lm2_m
443 angpart=dcos(md*phi+phase(l,m))
446 bpf(iph,ith)=bpf(iph,ith)+bamp(l,m)*angpart(iph)*cp_l_m(ith)
455 subroutine fft(a,b,ntot,n,nspan,isn)
499 double precision :: a(:),b(:)
503 dimension nfac(11),np(209)
505 dimension at(23),ck(23),bt(23),sk(23)
506 integer :: i,ii,maxp,maxf,n,inc,isn,nt,ntot,ks,nspan,kspan,nn,jc,jf,m
507 integer :: k,j,jj,nfac,kt,np,kk,k1,k2,k3,k4,kspnn
508 double precision :: c72,s72,s120,rad,radf,sd,cd,ak,bk,c1
509 double precision :: s1,aj,bj,akp,ajp,ajm,akm,bkp,bkm,bjp,bjm,aa
510 double precision :: bb,sk,ck,at,bt,s3,c3,s2,c2
524 c72=0.30901699437494742d0
525 s72=0.95105651629515357d0
526 s120=0.86602540378443865d0
527 rad=6.2831853071796d0
528 if(isn .ge. 0)
go to 10
538 radf=rad*dble(jc)*0.5d0
548 20
if(k-(k/16)*16 .eq. 0)
go to 15
555 30
if(mod(k,jj) .eq. 0)
go to 25
558 if(jj .le. k)
go to 30
559 if(k .gt. 4)
go to 40
564 40
if(k-(k/4)*4 .ne. 0)
go to 50
570 60
if(mod(k,j) .ne. 0)
go to 70
575 if(j .le. k)
go to 60
576 80
if(kt .eq. 0)
go to 100
581 if(j .ne. 0)
go to 90
583 100 sd=radf/dble(kspan)
588 if(nfac(i) .ne. 2)
go to 400
600 if(kk .le. nn)
go to 210
602 if(kk .le. jc)
go to 210
603 if(kk .gt. kspan)
go to 800
614 if(kk .lt. nt)
go to 230
618 if(kk .gt. k2)
go to 230
621 c1=2.d0-(ak**2+s1**2)
625 if(kk .lt. k2)
go to 230
628 if(kk .le. jc+jc)
go to 220
641 aj=(a(k1)-a(k2))*s120
642 bj=(b(k1)-b(k2))*s120
648 if(kk .lt. nn)
go to 320
650 if(kk .le. kspan)
go to 320
653 400
if(nfac(i) .ne. 4)
go to 600
673 if(isn .lt. 0)
go to 450
678 if(s1 .eq. 0.d0)
go to 460
679 430 a(k1)=akp*c1-bkp*s1
686 if(kk .le. nt)
go to 420
687 440 c2=c1-(cd*c1+sd*s1)
689 c1=2.d0-(c2**2+s1**2)
697 if(kk .le. kspan)
go to 420
699 if(kk .le. jc)
go to 410
700 if(kspan .eq. jc)
go to 800
706 if(s1 .ne. 0)
go to 430
714 if(kk .le. nt)
go to 420
752 if(kk .lt. nn)
go to 520
754 if(kk .le. kspan)
go to 520
760 if(k .eq. 3)
go to 320
761 if(k .eq. 5)
go to 510
762 if(k .eq. jf)
go to 640
767 if(jf .gt. maxf)
go to 998
771 630 ck(j)=ck(k)*c1+sk(k)*s1
772 sk(j)=ck(k)*s1-sk(k)*c1
777 if(j .lt. k)
go to 630
796 if(k1 .lt. k2)
go to 650
817 if(jj .gt. jf) jj=jj-jf
818 if(k .lt. jf)
go to 670
825 if(j .lt. k)
go to 660
827 if(kk .le. nn)
go to 640
829 if(kk .le. kspan)
go to 640
831 700
if(i .eq. m)
go to 800
842 if(kk .le. nt)
go to 730
847 if(kk .le. kspnn)
go to 730
850 c1=2.d0-(c2**2+s1**2)
854 if(kk .le. kspan)
go to 720
856 if(kk .le. jc+jc)
go to 710
861 if(kt .eq. 0)
go to 890
866 810 np(j+1)=np(j)/nfac(j)
867 np(k)=np(k+1)*nfac(j)
870 if(j .lt. k)
go to 810
876 if(n .ne. ntot)
go to 850
886 if(k2 .lt. ks)
go to 820
890 if(k2 .gt. np(j))
go to 830
892 840
if(kk .lt. k2)
go to 820
895 if(k2 .lt. ks)
go to 840
896 if(kk .lt. ks)
go to 830
909 if(kk .lt. k)
go to 860
912 if(kk .lt. nt)
go to 850
915 if(k2 .lt. ks)
go to 850
919 if(k2 .gt. np(j))
go to 870
921 880
if(kk .lt. k2)
go to 850
924 if(k2 .lt. ks)
go to 880
925 if(kk .lt. ks)
go to 870
927 890
if(2*kt+1 .ge. m)
return
932 900 nfac(j)=nfac(j)*nfac(j+1)
934 if(j .ne. kt)
go to 900
937 if(nn .gt. maxp)
go to 998
946 if(jj .ge. k2)
go to 902
952 if(j .le. nn)
go to 904
959 if(kk .ne. j)
go to 910
963 if(kk .lt. 0)
go to 914
964 if(kk .ne. j)
go to 910
966 if(j .ne. nn)
go to 914
971 if(np(j) .lt. 0)
go to 924
974 if(jj .gt. maxf) kspan=maxf
984 if(k1 .ne. kk)
go to 928
992 if(k1 .ne. kk)
go to 936
994 if(k .ne. j)
go to 932
1001 if(k1 .ne. kk)
go to 940
1002 if(jj .ne. 0)
go to 926
1003 if(j .ne. 1)
go to 924
1007 if(nt .ge. 0)
go to 924
1013 999
format(44h0array bounds exceeded within
subroutine fft)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
module mod_pfss.t – potential field source surface model PURPOSE : to extrapolate global potential ma...
subroutine fft(a, b, ntot, n, nspan, isn)
subroutine coef(b_r0, xm, ym, miu, mius, cfwm)
subroutine, public pfss(ixIL, ixOL, Bpf, x)
subroutine inv_sph_transform(Bt, phi, miu, mius, nphi, ntheta, Bpf, qlmax)
double precision, public r_0
subroutine cfweights(ym, miu, cfwm)
subroutine, public harm_coef(mapname)
double precision, public r_s