MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_pfss.t
Go to the documentation of this file.
1 !> module mod_pfss.t -- potential field source surface model
2 !> PURPOSE : to extrapolate global potential magnetic field of the sun from
3 !> synoptic magnetograms
4 !> 2013.11.04 Developed by S. Moschou and C. Xia
5 !> 2014.04.01 Allow to change source surface (C. Xia)
6 !> PRECONDITIONS:
7 !> 1. 3D spherical coordinates
8 !> 2. A synoptic magnetogram in a binary file contains nphi, ntheta,
9 !> theta(ntheta), phi(nphi), B_r(nphi,ntheta) succesively.
10 !> 3. nphi, ntheta are long integers and other arrays are double precision.
11 !> theta contains decreasing radians with increasing indice (Pi to 0)
12 !> phi contains increasing radians with increasing indice (0 to 2*Pi)
13 !> USAGE:
14 !> example for a magnetogram with name 'mdicr2020.dat':
15 !>
16 !> subroutine initglobaldata_usr
17 !> ...
18 !> R_0=1.d0 ! dimensionless Solar radius (default 1.0)
19 !> R_s=2.5d0 ! dimensionless radius of source surface (default 2.5)
20 !> lmax=120 ! use a fixed value instead of the value determined by the
21 !> resolution of input magnetogram
22 !> trunc=.true. ! use less spherical harmonics at larger distances
23 !> call harm_coef('mdicr2020.dat')
24 !> end subroutine initglobaldata_usr
25 !>
26 !> subroutine initonegrid_usr(ixI^L,ixO^L,w,x)
27 !> ...
28 !> double precision :: bpf(ixI^S,1:ndir)
29 !> ...
30 !> call pfss(ixI^L,ixO^L,bpf,x)
31 !> w(ix^S,mag(:))=bpf(ix^S,:)
32 !>
33 !> end subroutine initonegrid_usr
34 module mod_pfss
35  implicit none
36  private
37 
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(:)
42  integer, public :: lmax=0
43  logical, public :: trunc=.false.
44 
45  public :: harm_coef
46 {^ifthreed
47  public :: pfss
48 }
49 
50 contains
51 
52  subroutine harm_coef(mapname)
54  use mod_comm_lib, only: mpistop
55 
56  double precision, allocatable :: b_r0(:,:)
57  double precision, allocatable :: theta(:),phi(:),cfwm(:)
58  double precision :: rsl,xrl,dxr
59  integer :: xm,ym,l,m,amode,file_handle,il,ir,nlarr,nsh
60  integer, dimension(MPI_STATUS_SIZE) :: statuss
61  character(len=*) :: mapname
62  character(len=80) :: fharmcoef
63  logical :: aexist
64 
65  fharmcoef=mapname//'.coef'
66  inquire(file=fharmcoef, exist=aexist)
67  if(aexist) then
68  if(mype==0) then
69  call mpi_file_open(mpi_comm_self,fharmcoef,mpi_mode_rdonly, &
70  mpi_info_null,file_handle,ierrmpi)
71  call mpi_file_read(file_handle,lmax,1,mpi_integer,statuss,ierrmpi)
72  allocate(flm(0:lmax,0:lmax))
73  call mpi_file_read(file_handle,flm,(lmax+1)*(lmax+1),&
74  mpi_double_complex,statuss,ierrmpi)
75  call mpi_file_close(file_handle,ierrmpi)
76  end if
77  call mpi_barrier(icomm,ierrmpi)
78  if(npe>0) call mpi_bcast(lmax,1,mpi_integer,0,icomm,ierrmpi)
79  if(mype/=0) allocate(flm(0:lmax,0:lmax))
80  call mpi_barrier(icomm,ierrmpi)
81  if(npe>0) call mpi_bcast(flm,(lmax+1)*(lmax+1),mpi_double_complex,0,icomm,&
82  ierrmpi)
83  else
84  if(mype==0) then
85  inquire(file=mapname,exist=aexist)
86  if(.not. aexist) then
87  if(mype==0) write(*,'(2a)') "can not find file:",mapname
88  call mpistop("no input magnetogram found")
89  end if
90  call mpi_file_open(mpi_comm_self,mapname,mpi_mode_rdonly,mpi_info_null,&
91  file_handle,ierrmpi)
92  call mpi_file_read(file_handle,xm,1,mpi_integer,statuss,ierrmpi)
93  call mpi_file_read(file_handle,ym,1,mpi_integer,statuss,ierrmpi)
94  if(lmax==0) lmax=min(2*ym/3,xm/3)
95 
96  allocate(b_r0(xm,ym))
97  allocate(theta(ym))
98  allocate(phi(xm))
99  call mpi_file_read(file_handle,theta,ym,mpi_double_precision,&
100  statuss,ierrmpi)
101  call mpi_file_read(file_handle,phi,xm,mpi_double_precision,&
102  statuss,ierrmpi)
103  call mpi_file_read(file_handle,b_r0,xm*ym,mpi_double_precision,&
104  statuss,ierrmpi)
105  call mpi_file_close(file_handle,ierrmpi)
106  print*,'nphi,ntheta',xm,ym
107  print*,'theta range:',minval(theta),maxval(theta)
108  print*,'phi range:',minval(phi),maxval(phi)
109  print*,'Brmax,Brmin',maxval(b_r0),minval(b_r0)
110  allocate(cfwm(ym))
111  call cfweights(ym,dcos(theta),cfwm)
112  allocate(flm(0:lmax,0:lmax))
113  call coef(b_r0,xm,ym,dcos(theta),dsin(theta),cfwm)
114  deallocate(b_r0)
115  deallocate(theta)
116  deallocate(phi)
117  amode=ior(mpi_mode_create,mpi_mode_wronly)
118  call mpi_file_open(mpi_comm_self,fharmcoef,amode, &
119  mpi_info_null,file_handle,ierrmpi)
120  call mpi_file_write(file_handle,lmax,1,mpi_integer,statuss,ierrmpi)
121  call mpi_file_write(file_handle,flm,(lmax+1)*(lmax+1),&
122  mpi_double_complex,statuss,ierrmpi)
123  call mpi_file_close(file_handle,ierrmpi)
124  endif
125  call mpi_barrier(icomm,ierrmpi)
126  if(npe>1) call mpi_bcast(lmax,1,mpi_integer,0,icomm,ierrmpi)
127  if(mype/=0) allocate(flm(0:lmax,0:lmax))
128  call mpi_barrier(icomm,ierrmpi)
129  if(npe>1) call mpi_bcast(flm,(lmax+1)*(lmax+1),mpi_double_complex,0,&
130  icomm,ierrmpi)
131  end if
132  if(mype==0) print*,'lmax=',lmax,'trunc=',trunc
133  nlarr=501
134  allocate(lmaxarray(nlarr))
135  allocate(xrg(nlarr))
136  lmaxarray=lmax
137  if(trunc) then
138  dxr=(r_s-r_0)/dble(nlarr-1)
139  do ir=1,nlarr
140  xrg(ir)=dxr*dble(ir-1)+r_0
141  do il=0,lmax
142  xrl=xrg(ir)**il
143  if(xrl > 1.d6) then
144  lmaxarray(ir)=il
145  exit
146  end if
147  end do
148  end do
149  endif
150  ! calculate global Alm Blm Rlm
151  allocate(alm(0:lmax,0:lmax))
152  allocate(blm(0:lmax,0:lmax))
153  allocate(rlm(0:lmax,0:lmax))
154  alm=(0.d0,0.d0)
155  blm=(0.d0,0.d0)
156  do l=0,lmax
157  do m=0,l
158  rsl=r_s**(-(2*l+1))
159  rlm(l,m)=dsqrt(dble(l**2-m**2)/dble(4*l**2-1))
160  blm(l,m)=-flm(l,m)/(1.d0+dble(l)+dble(l)*rsl)
161  alm(l,m)=-rsl*blm(l,m)
162  end do
163  end do
164 
165  end subroutine harm_coef
166 
167  subroutine cfweights(ym,miu,cfwm)
169 
170  integer, intent(in) :: ym
171  double precision, intent(in) :: miu(ym)
172  double precision, intent(out) :: cfwm(ym)
173 
174  double precision,dimension(ym) :: Pl,Pm2,Pm1,Pprime,sintheta
175  double precision :: lr
176  integer :: l
177 
178  sintheta=dsqrt(1.d0-miu**2)
179 
180  pm2=1.d0
181  pm1=miu
182 
183  do l=2,ym-1
184  lr=1.d0/dble(l)
185  pl=(2.d0-lr)*pm1*miu-(1.d0-lr)*pm2
186  pm2=pm1
187  pm1=pl
188  end do
189 
190  pprime=(dble(ym)*pl)/sintheta**2
191  cfwm=2.d0/(sintheta*pprime)**2
192  cfwm=cfwm*(2.d0*dpi)
193 
194  end subroutine cfweights
195 
196  subroutine coef(b_r0,xm,ym,miu,mius,cfwm)
198 
199  integer, intent(in) :: xm,ym
200  double precision, intent(in) :: b_r0(xm,ym),cfwm(ym),miu(ym),mius(ym)
201 
202  double complex :: Bm(0:xm-1,0:ym-1)
203  double precision,dimension(xm) :: fftmr,fftmi
204  double precision,dimension(0:lmax) :: N_mm
205  double precision,dimension(ym) :: P_lm1,P_lm2,old_Pmm,P_l
206  double precision :: mr,lr,c1,c2
207  integer :: l,m,i,j,stat
208 
209  bm=(0.d0,0.d0)
210  do i=1,ym
211  fftmr=b_r0(:,i)/dble(xm)
212  fftmi=0.d0
213  call fft(fftmr,fftmi,xm,xm,xm,-1)
214  bm(:,i-1)=(fftmr+(0.d0,1.d0)*fftmi)
215  end do
216  n_mm(0)=1.d0/dsqrt(4.d0*dpi)
217  do m=1,lmax
218  n_mm(m)=-n_mm(m-1)*dsqrt(1.d0+1.d0/dble(2*m))
219  end do
220  !first do m=0
221  p_lm2=n_mm(0)
222  p_lm1=p_lm2*miu*dsqrt(3.d0)
223  !set l=0 m=0 term
224  flm(0,0)=sum(bm(0,:)*p_lm2*cfwm)
225  !set l=1 m=0 term
226  flm(1,0)=sum(bm(0,:)*p_lm1*cfwm)
227  do l=2,lmax
228  lr=dble(l)
229  c1=dsqrt(4.d0-1.d0/lr**2)
230  c2=-(1.d0-1.d0/lr)*dsqrt((2.d0*lr+1.d0)/(2.d0*lr-3.d0))
231  p_l=c1*miu*p_lm1+c2*p_lm2
232  !set m=0 term for all other l's
233  flm(l,0)=sum(bm(0,:)*p_l*cfwm)
234  p_lm2=p_lm1
235  p_lm1=p_l
236  end do
237 
238  !since only l modes from 0 to lmax are used
239  bm=2.d0*bm
240 
241  !now the rest of the m's
242  old_pmm=n_mm(0)
243  do m=1,lmax
244  p_lm2=old_pmm*mius*n_mm(m)/n_mm(m-1)
245  p_lm1=p_lm2*miu*dsqrt(dble(2*m+3))
246  !ACCURATE UP TO HERE
247  old_pmm=p_lm2
248  !set l=m mode
249  flm(m,m)=sum(bm(m,:)*p_lm2*cfwm)
250  !set l=m+1 mode
251  if(m<lmax) flm(m+1,m)=sum(bm(m,:)*p_lm1*cfwm)
252  mr=dble(m)
253  do l=m+2,lmax
254  lr=dble(l)
255  c1=dsqrt((4.d0*lr**2-1.d0)/(lr**2-mr**2))
256  c2=-dsqrt(((2.d0*lr+1.d0)*((lr-1.d0)**2-mr**2))/((2.d0*lr-3.d0)*(lr**2-&
257  mr**2)))
258  p_l=c1*miu*p_lm1+c2*p_lm2
259  flm(l,m)=sum(bm(m,:)*p_l*cfwm)
260  p_lm2=p_lm1
261  p_lm1=p_l
262  end do
263  end do
264 
265  end subroutine coef
266 {^ifthreed
267  subroutine pfss(ixI^L,ixO^L,Bpf,x)
269 
270  integer, intent(in) :: ixi^l,ixo^l
271  double precision, intent(in) :: x(ixi^s,1:ndim)
272  double precision, intent(out) :: bpf(ixi^s,1:ndir)
273 
274  double complex :: bt(0:lmax,0:lmax,ixomin1:ixomax1)
275  double precision :: phase(ixi^s,1:ndir),bpfiv(ixomin3:ixomax3,ixomin2:ixomax2)
276  double precision :: miu(ixomin2:ixomax2),mius(ixomin2:ixomax2),xr
277  integer :: l,m,ix^d,j,l1,l2,ntheta,nphi,ir,qlmax
278 
279  bt=(0.d0,0.d0)
280  nphi=ixomax3-ixomin3+1
281  ntheta=ixomax2-ixomin2+1
282  miu(ixomin2:ixomax2)=dcos(x(ixomin1,ixomax2:ixomin2:-1,ixomin3,2))
283  mius(ixomin2:ixomax2)=dsin(x(ixomin1,ixomax2:ixomin2:-1,ixomin3,2))
284  do ix1=ixomin1,ixomax1
285  xr=x(ix1,ixomin2,ixomin3,1)
286  if(trunc) then
287  do ir=1,size(lmaxarray)
288  if(xrg(ir)>=xr) exit
289  end do
290  if(ir>size(lmaxarray)) ir=size(lmaxarray)
291  qlmax=lmaxarray(ir)
292  else
293  qlmax=lmax
294  endif
295  !Calculate Br
296  do l=0,lmax
297  do m=0,l
298  bt(l,m,ix1)=alm(l,m)*dble(l)*xr**(l-1)-blm(l,m)*dble(l+1)*xr**(-l-2)
299  end do
300  enddo
301  call inv_sph_transform(bt(:,:,ix1),x(ixomin1,ixomin2,&
302  ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
303  do ix3=ixomin3,ixomax3
304  do ix2=ixomin2,ixomax2
305  bpf(ix1,ix2,ix3,1)=bpfiv(ix3,ixomax2-ix2+ixomin2)
306  enddo
307  enddo
308  !Calculate Btheta
309  do l=0,lmax
310  do m=0,l
311  if (l==0) then
312  bt(l,m,ix1)=-rlm(l+1,m)*dble(l+2)*&
313  (alm(l+1,m)*xr**l+blm(l+1,m)*xr**(-l-3))
314  else if (l>=1 .and. l<=lmax-1) then
315  bt(l,m,ix1)=rlm(l,m)*&
316  dble(l-1)*(alm(l-1,m)*xr**(l-2)+blm(l-1,m)*&
317  xr**(-l-1))-rlm(l+1,m)*dble(l+2)*&
318  (alm(l+1,m)*xr**l+blm(l+1,m)*xr**(-l-3))
319  else
320  bt(l,m,ix1)=rlm(l,m)*&
321  dble(l-1)*(alm(l-1,m)*xr**(l-2)+blm(l-1,m)*xr**(-l-1))
322  end if
323  end do
324  enddo
325  call inv_sph_transform(bt(:,:,ix1),x(ixomin1,ixomin2,&
326  ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
327  do ix3=ixomin3,ixomax3
328  do ix2=ixomin2,ixomax2
329  bpf(ix1,ix2,ix3,2)=bpfiv(ix3,ixomax2-ix2+ixomin2)/mius(&
330  ixomax2-ix2+ixomin2)
331  enddo
332  enddo
333 
334  !Calculate Bphi
335  do l=0,lmax
336  do m=0,l
337  bt(l,m,ix1)=(0.d0,1.d0)*m*(alm(l,m)*xr**(l-1)+blm(l,m)*xr**(-l-2))
338  end do
339  enddo
340  call inv_sph_transform(bt(:,:,ix1),x(ixomin1,ixomin2,&
341  ixomin3:ixomax3,3),miu,mius,nphi,ntheta,bpfiv,qlmax)
342  do ix3=ixomin3,ixomax3
343  do ix2=ixomin2,ixomax2
344  bpf(ix1,ix2,ix3,3)=bpfiv(ix3,ixomax2-ix2+ixomin2)/mius(&
345  ixomax2-ix2+ixomin2)
346  enddo
347  enddo
348  enddo
349  !Scalar Potential
350  ! Potlc(ix^D)=Alm(l,m)*x(ix^D,1)**l+Blm(l,m)*x(ix^D,1)**(-l-1)
351 
352  !do ix3=ixOmin3,ixOmax3
353  ! do ix2=ixOmin2,ixOmax2
354  ! print*,x(ix1,ixOmax2-ix2+ixOmin2,ix3,2)
355  ! print*,'miu==',miu(ixOmax2-ix2+ixOmin2)
356  ! enddo
357  !enddo
358 
359  end subroutine pfss
360 
361  subroutine inv_sph_transform(Bt,phi,miu,mius,nphi,ntheta,Bpf,qlmax)
363 
364  integer, intent(in) :: nphi,ntheta,qlmax
365  double complex, intent(in) :: Bt(0:lmax,0:lmax)
366  double precision, intent(in) :: phi(nphi),miu(ntheta),mius(ntheta)
367  double precision, intent(out) :: Bpf(nphi,ntheta)
368 
369  double precision,dimension(0:lmax,0:lmax) :: cp,phase,Bamp
370  double precision,dimension(ntheta) :: cp_1_0,cp_l_0,cp_lm1_0,cp_lm2_0,cp_m_m
371  double precision,dimension(ntheta) :: cp_1_m,cp_l_m,cp_lm1_m,cp_lm2_m,cp_mp1_m
372  double precision :: angpart(nphi)
373  double precision :: ld,md,c1,c2,cp_0_0
374  integer :: l,m,iph,ith
375 
376  bamp=abs(bt)
377 
378  phase=atan2(dimag(bt),dble(bt))
379 
380  bpf=0.d0
381  !take care of modes where m=0
382  cp_0_0=dsqrt(1.d0/(4.d0*dpi))
383  !start with l=m=0 mode
384  bpf=bpf+bamp(0,0)*dcos(phase(0,0))*cp_0_0
385 
386  !proceed with l=1 m=0 mode
387  cp_1_0=dsqrt(3.d0)*miu*cp_0_0
388  do iph=1,nphi
389  bpf(iph,:)=bpf(iph,:)+bamp(1,0)*dcos(phase(1,0))*cp_1_0
390  enddo
391 
392  !proceed with l modes for which m=0
393  cp_lm1_0=cp_0_0
394  cp_l_0=cp_1_0
395  do l=2,qlmax
396  ld=dble(l)
397  cp_lm2_0=cp_lm1_0
398  cp_lm1_0=cp_l_0
399  c1=dsqrt(4.d0*ld**2-1.d0)/ld
400  c2=dsqrt((2.d0*ld+1.d0)/(2.d0*ld-3.d0))*(ld-1.d0)/ld
401  cp_l_0=c1*miu*cp_lm1_0-c2*cp_lm2_0
402  do iph=1,nphi
403  bpf(iph,:)=bpf(iph,:)+bamp(l,0)*dcos(phase(l,0))*cp_l_0
404  enddo
405  enddo
406 
407  !loop through m's for m>0 and then loop through l's for each m
408  cp_m_m=cp_0_0
409  do m=1,qlmax
410  md=dble(m)
411  !first do l=m modes
412  cp_m_m=-dsqrt(1.d0+1.d0/(2.d0*md))*mius*cp_m_m
413  do iph=1,nphi
414  angpart(iph)=dcos(md*phi(iph)+phase(m,m))
415  end do
416  do ith=1,ntheta
417  do iph=1,nphi
418  bpf(iph,ith)=bpf(iph,ith)+bamp(m,m)*angpart(iph)*cp_m_m(ith)
419  enddo
420  enddo
421 
422  !proceed with l=m+1 modes
423  if(qlmax>=m+1) then
424  cp_mp1_m=dsqrt(2.d0*md+3.d0)*miu*cp_m_m
425  angpart=dcos(md*phi+phase(m+1,m))
426  do ith=1,ntheta
427  do iph=1,nphi
428  bpf(iph,ith)=bpf(iph,ith)+bamp(m+1,m)*angpart(iph)*cp_mp1_m(ith)
429  enddo
430  enddo
431  endif
432 
433  !finish with the rest l
434  if(qlmax>=m+2) then
435  cp_lm1_m=cp_m_m
436  cp_l_m=cp_mp1_m
437  do l=m+2,qlmax
438  ld=dble(l)
439  cp_lm2_m=cp_lm1_m
440  cp_lm1_m=cp_l_m
441  c1=dsqrt((4.d0*ld**2-1.d0)/(ld**2-md**2))
442  c2=dsqrt((2.d0*ld+1.d0)*((ld-1.d0)**2-md**2)/(2.d0*ld-3.d0)/(ld**2-md**2))
443  cp_l_m=c1*miu*cp_lm1_m-c2*cp_lm2_m
444  angpart=dcos(md*phi+phase(l,m))
445  do ith=1,ntheta
446  do iph=1,nphi
447  bpf(iph,ith)=bpf(iph,ith)+bamp(l,m)*angpart(iph)*cp_l_m(ith)
448  enddo
449  enddo
450  enddo
451  endif
452  enddo
453 
454  end subroutine inv_sph_transform
455 }
456  subroutine fft(a,b,ntot,n,nspan,isn)
457  ! multivariate complex fourier transform, computed in place
458  ! using mixed-radix fast fourier transform algorithm.
459  ! by r. c. singleton, stanford research institute, sept. 1968
460  ! arrays a and b originally hold the real and imaginary
461  ! components of the data, and return the real and
462  ! imaginary components of the resulting fourier coefficients.
463  ! multivariate data is indexed according to the fortran
464  ! array element successor function, without limit
465  ! on the number of implied multiple subscripts.
466  ! the subroutine is called once for each variate.
467  ! the calls for a multivariate transform may be in any order.
468  ! ntot is the total number of complex data values.
469  ! n is the dimension of the current variable.
470  ! nspan/n is the spacing of consecutive data values
471  ! while indexing the current variable.
472  ! the sign of isn determines the sign of the complex
473  ! exponential, and the magnitude of isn is normally one.
474  ! a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
475  ! is computed by
476  ! call fft(a,b,n1*n2*n3,n1,n1,1)
477  ! call fft(a,b,n1*n2*n3,n2,n1*n2,1)
478  ! call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
479  ! for a single-variate transform,
480  ! ntot = n = nspan = (number of complex data values), e.g.
481  ! call fft(a,b,n,n,n,1)
482  ! the data can alternatively be stored in a single complex array c
483  ! in standard fortran fashion, i.e. alternating real and imaginary
484  ! parts. then with most fortran compilers, the complex array c can
485  ! be equivalenced to a real array a, the magnitude of isn changed
486  ! to two to give correct indexing increment, and a(1) and a(2) used
487  ! to pass the initial addresses for the sequences of real and
488  ! imaginary values, e.g.
489  ! complex c(ntot)
490  ! real a(2*ntot)
491  ! equivalence (c(1),a(1))
492  ! call fft(a(1),a(2),ntot,n,nspan,2)
493  ! arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
494  ! are used for temporary storage. if the available storage
495  ! is insufficient, the program is terminated by a stop.
496  ! maxf must be .ge. the maximum prime factor of n.
497  ! maxp must be .gt. the number of prime factors of n.
498  ! in addition, if the square-free portion k of n has two or
499  ! more prime factors, then maxp must be .ge. k-1.
500  double precision :: a(:),b(:)
501  ! array storage in nfac for a maximum of 15 prime factors of n.
502  ! if n has more than one square-free factor, the product of the
503  ! square-free factors must be .le. 210
504  dimension nfac(11),np(209)
505  ! array storage for maximum prime factor of 23
506  dimension at(23),ck(23),bt(23),sk(23)
507  integer :: i,ii,maxp,maxf,n,inc,isn,nt,ntot,ks,nspan,kspan,nn,jc,jf,m
508  integer :: k,j,jj,nfac,kt,np,kk,k1,k2,k3,k4,kspnn
509  double precision :: c72,s72,s120,rad,radf,sd,cd,ak,bk,c1
510  double precision :: s1,aj,bj,akp,ajp,ajm,akm,bkp,bkm,bjp,bjm,aa
511  double precision :: bb,sk,ck,at,bt,s3,c3,s2,c2
512  equivalence(i,ii)
513 
514  ! the following two constants should agree with the array dimensions.
515  maxp=209
516  ! Date: Wed, 9 Aug 1995 09:38:49 -0400
517  ! From: ldm@apollo.numis.nwu.edu
518  maxf=23
519  s3=0.d0
520  s2=0.d0
521  c3=0.d0
522  c2=0.d0
523  if(n .lt. 2) return
524  inc=isn
525  c72=0.30901699437494742d0
526  s72=0.95105651629515357d0
527  s120=0.86602540378443865d0
528  rad=6.2831853071796d0
529  if(isn .ge. 0) go to 10
530  s72=-s72
531  s120=-s120
532  rad=-rad
533  inc=-inc
534  10 nt=inc*ntot
535  ks=inc*nspan
536  kspan=ks
537  nn=nt-inc
538  jc=ks/n
539  radf=rad*dble(jc)*0.5d0
540  i=0
541  jf=0
542  ! determine the factors of n
543  m=0
544  k=n
545  go to 20
546  15 m=m+1
547  nfac(m)=4
548  k=k/16
549  20 if(k-(k/16)*16 .eq. 0) go to 15
550  j=3
551  jj=9
552  go to 30
553  25 m=m+1
554  nfac(m)=j
555  k=k/jj
556  30 if(mod(k,jj) .eq. 0) go to 25
557  j=j+2
558  jj=j**2
559  if(jj .le. k) go to 30
560  if(k .gt. 4) go to 40
561  kt=m
562  nfac(m+1)=k
563  if(k .ne. 1) m=m+1
564  go to 80
565  40 if(k-(k/4)*4 .ne. 0) go to 50
566  m=m+1
567  nfac(m)=2
568  k=k/4
569  50 kt=m
570  j=2
571  60 if(mod(k,j) .ne. 0) go to 70
572  m=m+1
573  nfac(m)=j
574  k=k/j
575  70 j=((j+1)/2)*2+1
576  if(j .le. k) go to 60
577  80 if(kt .eq. 0) go to 100
578  j=kt
579  90 m=m+1
580  nfac(m)=nfac(j)
581  j=j-1
582  if(j .ne. 0) go to 90
583  ! compute fourier transform
584  100 sd=radf/dble(kspan)
585  cd=2.d0*dsin(sd)**2
586  sd=dsin(sd+sd)
587  kk=1
588  i=i+1
589  if(nfac(i) .ne. 2) go to 400
590  ! transform for factor of 2 (including rotation factor)
591  kspan=kspan/2
592  k1=kspan+2
593  210 k2=kk+kspan
594  ak=a(k2)
595  bk=b(k2)
596  a(k2)=a(kk)-ak
597  b(k2)=b(kk)-bk
598  a(kk)=a(kk)+ak
599  b(kk)=b(kk)+bk
600  kk=k2+kspan
601  if(kk .le. nn) go to 210
602  kk=kk-nn
603  if(kk .le. jc) go to 210
604  if(kk .gt. kspan) go to 800
605  220 c1=1.d0-cd
606  s1=sd
607  230 k2=kk+kspan
608  ak=a(kk)-a(k2)
609  bk=b(kk)-b(k2)
610  a(kk)=a(kk)+a(k2)
611  b(kk)=b(kk)+b(k2)
612  a(k2)=c1*ak-s1*bk
613  b(k2)=s1*ak+c1*bk
614  kk=k2+kspan
615  if(kk .lt. nt) go to 230
616  k2=kk-nt
617  c1=-c1
618  kk=k1-k2
619  if(kk .gt. k2) go to 230
620  ak=c1-(cd*c1+sd*s1)
621  s1=(sd*c1-cd*s1)+s1
622  c1=2.d0-(ak**2+s1**2)
623  s1=c1*s1
624  c1=c1*ak
625  kk=kk+jc
626  if(kk .lt. k2) go to 230
627  k1=k1+inc+inc
628  kk=(k1-kspan)/2+jc
629  if(kk .le. jc+jc) go to 220
630  go to 100
631  ! transform for factor of 3 (optional code)
632  320 k1=kk+kspan
633  k2=k1+kspan
634  ak=a(kk)
635  bk=b(kk)
636  aj=a(k1)+a(k2)
637  bj=b(k1)+b(k2)
638  a(kk)=ak+aj
639  b(kk)=bk+bj
640  ak=-0.5d0*aj+ak
641  bk=-0.5d0*bj+bk
642  aj=(a(k1)-a(k2))*s120
643  bj=(b(k1)-b(k2))*s120
644  a(k1)=ak-bj
645  b(k1)=bk+aj
646  a(k2)=ak+bj
647  b(k2)=bk-aj
648  kk=k2+kspan
649  if(kk .lt. nn) go to 320
650  kk=kk-nn
651  if(kk .le. kspan) go to 320
652  go to 700
653  ! transform for factor of 4
654  400 if(nfac(i) .ne. 4) go to 600
655  kspnn=kspan
656  kspan=kspan/4
657  410 c1=1.d0
658  s1=0.d0
659  420 k1=kk+kspan
660  k2=k1+kspan
661  k3=k2+kspan
662  akp=a(kk)+a(k2)
663  akm=a(kk)-a(k2)
664  ajp=a(k1)+a(k3)
665  ajm=a(k1)-a(k3)
666  a(kk)=akp+ajp
667  ajp=akp-ajp
668  bkp=b(kk)+b(k2)
669  bkm=b(kk)-b(k2)
670  bjp=b(k1)+b(k3)
671  bjm=b(k1)-b(k3)
672  b(kk)=bkp+bjp
673  bjp=bkp-bjp
674  if(isn .lt. 0) go to 450
675  akp=akm-bjm
676  akm=akm+bjm
677  bkp=bkm+ajm
678  bkm=bkm-ajm
679  if(s1 .eq. 0.d0) go to 460
680  430 a(k1)=akp*c1-bkp*s1
681  b(k1)=akp*s1+bkp*c1
682  a(k2)=ajp*c2-bjp*s2
683  b(k2)=ajp*s2+bjp*c2
684  a(k3)=akm*c3-bkm*s3
685  b(k3)=akm*s3+bkm*c3
686  kk=k3+kspan
687  if(kk .le. nt) go to 420
688  440 c2=c1-(cd*c1+sd*s1)
689  s1=(sd*c1-cd*s1)+s1
690  c1=2.d0-(c2**2+s1**2)
691  s1=c1*s1
692  c1=c1*c2
693  c2=c1**2-s1**2
694  s2=2.d0*c1*s1
695  c3=c2*c1-s2*s1
696  s3=c2*s1+s2*c1
697  kk=kk-nt+jc
698  if(kk .le. kspan) go to 420
699  kk=kk-kspan+inc
700  if(kk .le. jc) go to 410
701  if(kspan .eq. jc) go to 800
702  go to 100
703  450 akp=akm+bjm
704  akm=akm-bjm
705  bkp=bkm-ajm
706  bkm=bkm+ajm
707  if(s1 .ne. 0) go to 430
708  460 a(k1)=akp
709  b(k1)=bkp
710  a(k2)=ajp
711  b(k2)=bjp
712  a(k3)=akm
713  b(k3)=bkm
714  kk=k3+kspan
715  if(kk .le. nt) go to 420
716  go to 440
717  ! transform for factor of 5 (optional code)
718  510 c2=c72**2-s72**2
719  s2=2.d0*c72*s72
720  520 k1=kk+kspan
721  k2=k1+kspan
722  k3=k2+kspan
723  k4=k3+kspan
724  akp=a(k1)+a(k4)
725  akm=a(k1)-a(k4)
726  bkp=b(k1)+b(k4)
727  bkm=b(k1)-b(k4)
728  ajp=a(k2)+a(k3)
729  ajm=a(k2)-a(k3)
730  bjp=b(k2)+b(k3)
731  bjm=b(k2)-b(k3)
732  aa=a(kk)
733  bb=b(kk)
734  a(kk)=aa+akp+ajp
735  b(kk)=bb+bkp+bjp
736  ak=akp*c72+ajp*c2+aa
737  bk=bkp*c72+bjp*c2+bb
738  aj=akm*s72+ajm*s2
739  bj=bkm*s72+bjm*s2
740  a(k1)=ak-bj
741  a(k4)=ak+bj
742  b(k1)=bk+aj
743  b(k4)=bk-aj
744  ak=akp*c2+ajp*c72+aa
745  bk=bkp*c2+bjp*c72+bb
746  aj=akm*s2-ajm*s72
747  bj=bkm*s2-bjm*s72
748  a(k2)=ak-bj
749  a(k3)=ak+bj
750  b(k2)=bk+aj
751  b(k3)=bk-aj
752  kk=k4+kspan
753  if(kk .lt. nn) go to 520
754  kk=kk-nn
755  if(kk .le. kspan) go to 520
756  go to 700
757  ! transform for odd factors
758  600 k=nfac(i)
759  kspnn=kspan
760  kspan=kspan/k
761  if(k .eq. 3) go to 320
762  if(k .eq. 5) go to 510
763  if(k .eq. jf) go to 640
764  jf=k
765  s1=rad/dble(k)
766  c1=dcos(s1)
767  s1=dsin(s1)
768  if(jf .gt. maxf) go to 998
769  ck(jf)=1.d0
770  sk(jf)=0.d0
771  j=1
772  630 ck(j)=ck(k)*c1+sk(k)*s1
773  sk(j)=ck(k)*s1-sk(k)*c1
774  k=k-1
775  ck(k)=ck(j)
776  sk(k)=-sk(j)
777  j=j+1
778  if(j .lt. k) go to 630
779  640 k1=kk
780  k2=kk+kspnn
781  aa=a(kk)
782  bb=b(kk)
783  ak=aa
784  bk=bb
785  j=1
786  k1=k1+kspan
787  650 k2=k2-kspan
788  j=j+1
789  at(j)=a(k1)+a(k2)
790  ak=at(j)+ak
791  bt(j)=b(k1)+b(k2)
792  bk=bt(j)+bk
793  j=j+1
794  at(j)=a(k1)-a(k2)
795  bt(j)=b(k1)-b(k2)
796  k1=k1+kspan
797  if(k1 .lt. k2) go to 650
798  a(kk)=ak
799  b(kk)=bk
800  k1=kk
801  k2=kk+kspnn
802  j=1
803  660 k1=k1+kspan
804  k2=k2-kspan
805  jj=j
806  ak=aa
807  bk=bb
808  aj=0.d0
809  bj=0.d0
810  k=1
811  670 k=k+1
812  ak=at(k)*ck(jj)+ak
813  bk=bt(k)*ck(jj)+bk
814  k=k+1
815  aj=at(k)*sk(jj)+aj
816  bj=bt(k)*sk(jj)+bj
817  jj=jj+j
818  if(jj .gt. jf) jj=jj-jf
819  if(k .lt. jf) go to 670
820  k=jf-j
821  a(k1)=ak-bj
822  b(k1)=bk+aj
823  a(k2)=ak+bj
824  b(k2)=bk-aj
825  j=j+1
826  if(j .lt. k) go to 660
827  kk=kk+kspnn
828  if(kk .le. nn) go to 640
829  kk=kk-nn
830  if(kk .le. kspan) go to 640
831  ! multiply by rotation factor (except for factors of 2 and 4)
832  700 if(i .eq. m) go to 800
833  kk=jc+1
834  710 c2=1.d0-cd
835  s1=sd
836  720 c1=c2
837  s2=s1
838  kk=kk+kspan
839  730 ak=a(kk)
840  a(kk)=c2*ak-s2*b(kk)
841  b(kk)=s2*ak+c2*b(kk)
842  kk=kk+kspnn
843  if(kk .le. nt) go to 730
844  ak=s1*s2
845  s2=s1*c2+c1*s2
846  c2=c1*c2-ak
847  kk=kk-nt+kspan
848  if(kk .le. kspnn) go to 730
849  c2=c1-(cd*c1+sd*s1)
850  s1=s1+(sd*c1-cd*s1)
851  c1=2.d0-(c2**2+s1**2)
852  s1=c1*s1
853  c2=c1*c2
854  kk=kk-kspnn+jc
855  if(kk .le. kspan) go to 720
856  kk=kk-kspan+jc+inc
857  if(kk .le. jc+jc) go to 710
858  go to 100
859  ! permute the results to normal order---done in two stages
860  ! permutation for square factors of n
861  800 np(1)=ks
862  if(kt .eq. 0) go to 890
863  k=kt+kt+1
864  if(m .lt. k) k=k-1
865  j=1
866  np(k+1)=jc
867  810 np(j+1)=np(j)/nfac(j)
868  np(k)=np(k+1)*nfac(j)
869  j=j+1
870  k=k-1
871  if(j .lt. k) go to 810
872  k3=np(k+1)
873  kspan=np(2)
874  kk=jc+1
875  k2=kspan+1
876  j=1
877  if(n .ne. ntot) go to 850
878  ! permutation for single-variate transform (optional code)
879  820 ak=a(kk)
880  a(kk)=a(k2)
881  a(k2)=ak
882  bk=b(kk)
883  b(kk)=b(k2)
884  b(k2)=bk
885  kk=kk+inc
886  k2=kspan+k2
887  if(k2 .lt. ks) go to 820
888  830 k2=k2-np(j)
889  j=j+1
890  k2=np(j+1)+k2
891  if(k2 .gt. np(j)) go to 830
892  j=1
893  840 if(kk .lt. k2) go to 820
894  kk=kk+inc
895  k2=kspan+k2
896  if(k2 .lt. ks) go to 840
897  if(kk .lt. ks) go to 830
898  jc=k3
899  go to 890
900  ! permutation for multivariate transform
901  850 k=kk+jc
902  860 ak=a(kk)
903  a(kk)=a(k2)
904  a(k2)=ak
905  bk=b(kk)
906  b(kk)=b(k2)
907  b(k2)=bk
908  kk=kk+inc
909  k2=k2+inc
910  if(kk .lt. k) go to 860
911  kk=kk+ks-jc
912  k2=k2+ks-jc
913  if(kk .lt. nt) go to 850
914  k2=k2-nt+kspan
915  kk=kk-nt+jc
916  if(k2 .lt. ks) go to 850
917  870 k2=k2-np(j)
918  j=j+1
919  k2=np(j+1)+k2
920  if(k2 .gt. np(j)) go to 870
921  j=1
922  880 if(kk .lt. k2) go to 850
923  kk=kk+jc
924  k2=kspan+k2
925  if(k2 .lt. ks) go to 880
926  if(kk .lt. ks) go to 870
927  jc=k3
928  890 if(2*kt+1 .ge. m) return
929  kspnn=np(kt+1)
930  ! permutation for square-free factors of n
931  j=m-kt
932  nfac(j+1)=1
933  900 nfac(j)=nfac(j)*nfac(j+1)
934  j=j-1
935  if(j .ne. kt) go to 900
936  kt=kt+1
937  nn=nfac(kt)-1
938  if(nn .gt. maxp) go to 998
939  jj=0
940  j=0
941  go to 906
942  902 jj=jj-k2
943  k2=kk
944  k=k+1
945  kk=nfac(k)
946  904 jj=kk+jj
947  if(jj .ge. k2) go to 902
948  np(j)=jj
949  906 k2=nfac(kt)
950  k=kt+1
951  kk=nfac(k)
952  j=j+1
953  if(j .le. nn) go to 904
954  ! determine the permutation cycles of length greater than 1
955  j=0
956  go to 914
957  910 k=kk
958  kk=np(k)
959  np(k)=-kk
960  if(kk .ne. j) go to 910
961  k3=kk
962  914 j=j+1
963  kk=np(j)
964  if(kk .lt. 0) go to 914
965  if(kk .ne. j) go to 910
966  np(j)=-j
967  if(j .ne. nn) go to 914
968  maxf=inc*maxf
969  ! reorder a and b, following the permutation cycles
970  go to 950
971  924 j=j-1
972  if(np(j) .lt. 0) go to 924
973  jj=jc
974  926 kspan=jj
975  if(jj .gt. maxf) kspan=maxf
976  jj=jj-kspan
977  k=np(j)
978  kk=jc*k+ii+jj
979  k1=kk+kspan
980  k2=0
981  928 k2=k2+1
982  at(k2)=a(k1)
983  bt(k2)=b(k1)
984  k1=k1-inc
985  if(k1 .ne. kk) go to 928
986  932 k1=kk+kspan
987  k2=k1-jc*(k+np(k))
988  k=-np(k)
989  936 a(k1)=a(k2)
990  b(k1)=b(k2)
991  k1=k1-inc
992  k2=k2-inc
993  if(k1 .ne. kk) go to 936
994  kk=k2
995  if(k .ne. j) go to 932
996  k1=kk+kspan
997  k2=0
998  940 k2=k2+1
999  a(k1)=at(k2)
1000  b(k1)=bt(k2)
1001  k1=k1-inc
1002  if(k1 .ne. kk) go to 940
1003  if(jj .ne. 0) go to 926
1004  if(j .ne. 1) go to 924
1005  950 j=k3+1
1006  nt=nt-kspnn
1007  ii=nt-inc+1
1008  if(nt .ge. 0) go to 924
1009  return
1010  ! error finish, insufficient array storage
1011  998 isn=0
1012  print 999
1013  stop
1014  999 format(44h0array bounds exceeded within subroutine fft)
1015  end subroutine fft
1016 
1017 end module mod_pfss
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
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...
Definition: mod_pfss.t:34
subroutine fft(a, b, ntot, n, nspan, isn)
Definition: mod_pfss.t:457
subroutine coef(b_r0, xm, ym, miu, mius, cfwm)
Definition: mod_pfss.t:197
subroutine, public pfss(ixIL, ixOL, Bpf, x)
Definition: mod_pfss.t:268
subroutine inv_sph_transform(Bt, phi, miu, mius, nphi, ntheta, Bpf, qlmax)
Definition: mod_pfss.t:362
double precision, public r_0
Definition: mod_pfss.t:40
subroutine cfweights(ym, miu, cfwm)
Definition: mod_pfss.t:168
integer, public lmax
Definition: mod_pfss.t:42
subroutine, public harm_coef(mapname)
Definition: mod_pfss.t:53
logical, public trunc
Definition: mod_pfss.t:43
double precision, public r_s
Definition: mod_pfss.t:40