PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
dtorh.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
4 !------------------------------------------------------------------------------!
7 !------------------------------------------------------------------------------!
8 module dtorh
9 #include <PB3D_macros.h>
10  use num_vars, only: dp, pi, max_str_ln
11  use messages
12  use str_utilities
13 
14  implicit none
15  private
16  public dtorh1
17 
18  character(len=max_str_ln) :: err_msg
19 
20 contains
21  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
22  ! INPUT : !
23  ! Z ARGUMENT OF THE FUNCTIONS !
24  ! M ORDER OF THE FUNCTIONS !
25  ! NMAX MAXIMUM DEGREE OF THE FUNCTIONS : !
26  ! WE GET FUNCTIONS OF ALL THE ORDERS BELOW !
27  ! MIN(NEWN,NMAX). NEWN IS DEFINED BELOW . !
28  ! OUTPUT : !
29  ! *IF MODE IS EQUAL TO 0: !
30  ! PL(N) !
31  ! THESE VALUES ARE KEPT IN AN ARRAY !
32  ! QL(N) !
33  ! THESE VALUES ARE KEPT IN AN ARRAY !
34  ! !
35  ! NEWN MAXIMUM ORDER OF FUNCTIONS CALCULATED WHEN !
36  ! PL (NMAX+1) IS LARGER THAN 1/TINY !
37  ! (OVERFLOW LIMIT = 1/TINY, TINY IS DEFINED BELOW) !
38  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74  integer function dtorh1(z,m,nmax,pl,ql,newn,mode,ipre) result(ierr)
75  character(*), parameter :: rout_name = 'DTORH1'
76 
77  ! input / output
78  real(dp), intent(in) :: z
79  integer, intent(in) :: m
80  integer, intent(in) :: nmax
81  real(dp), intent(inout) :: pl(0:nmax)
82  real(dp), intent(inout) :: ql(0:nmax)
83  integer, intent(inout) :: newn
84  integer, intent(in), optional :: mode
85  integer, intent(in), optional :: ipre
86 
87  ! local variables
88  INTEGER nmaxp,mp,np,n,ical,i
89  INTEGER :: mode_loc, ipre_loc
90  REAL(dp) over,tinysq,qz,pisq,fl,cc,ar,gamma,fc,qdc1,qargu,&
91  &ARGU1,DFACQS,FCP,DD,D1,QM0,DFAC3,DFAC4,PL0
92  real(dp), allocatable :: pl_loc(:), ql_loc(:) ! local copy of pl and ql, possibly larger
93  !REAL(dp) DPPI ! Is never used?
94  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95  ! THE DIMENSION OF QLMM (INTERNAL ARRAY) MUST BE GREATER THAN M !
96  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97  REAL(dp) qlmm(0:1001),pr(2)
98  real(dp), PARAMETER :: eps=1.e-14_dp, tiny=1.e-290_dp
99 
100  ! initialize ierr
101  ierr = 0
102 
103  ! test
104  if (nmax.lt.0) then
105  ierr = 1
106  err_msg = 'Lowest possible degree is 0. Use recursive properties &
107  &of toroidal functions if you need lower.'
108  chckerr(err_msg)
109  end if
110 
111  ! set optional parameters
112  mode_loc = 0
113  if (present(mode)) mode_loc = mode
114  ipre_loc = 1
115  if (present(ipre)) ipre_loc = ipre
116 
117  over=1._dp/tiny
118  tinysq=dsqrt(tiny)
119  IF ((ipre_loc.NE.1).AND.(ipre_loc.NE.2)) THEN
120  ierr = 1
121  err_msg = 'IPRE MUST BE 1 OR 2'
122  chckerr(err_msg)
123  END IF
124 
125  pr(1)=.22_dp
126  pr(2)=.12_dp
127  nmaxp=nmax-1 ! decremented by one to get more logical ouptut
128  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129  ! EPS: REQUIRED ACCURACY FOR THE CONTINUED FRACTION !
130  ! (MODIFIED LENTZ) !
131  ! TINY: SMALL PARAMETER TO PREVENT OVERFLOWS IN THE CF !
132  ! (CLOSE TO THE UNDERFLOW LIMIT) !
133  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
134  IF (z.LE.1._dp) THEN
135  ierr = 1
136  err_msg = 'IMPROPER ARGUMENT. Z MUST BE GREATER THAN 1'
137  chckerr(err_msg)
138  END IF
139  qz=z
140  pisq=dsqrt(pi)
141  !DPPI=DSQRT(2._dp)/PISQ ! Is never used?
142  fl=m/2._dp
143  cc=abs(float(int(fl))-fl)
144  IF (cc.LT.0.4_dp) THEN
145  ar=1.
146  ELSE
147  ar=-1.
148  END IF
149  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150  ! WE CHOOSE EXPANSION OR CF FOR PL(0) DEPENDING ON THE VALUES !
151  ! OF Z,M AND MODE !
152  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153  ical=1
154  IF (m.NE.0) THEN
155  IF ((z/m).GT.pr(ipre_loc)) ical=2
156  IF (z.LT.5._dp) THEN
157  ical=1
158  ELSE IF (z.GT.20.d0) THEN
159  IF ((mode_loc.NE.2).AND.(ical.EQ.1)) ical=0
160  END IF
161  IF (ical.EQ.0) THEN
162  ierr =1
163  chckerr('YOU MUST CHOOSE MODE=2')
164  END IF
165  ELSE
166  IF (z.LT.5.d0) THEN
167  ical=1
168  ELSE
169  ical=2
170  END IF
171  END IF
172  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173  ! WE USE THE CODE IF NMAX IS GREATER THAN OR EQUAL TO 2 !
174  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175  IF(nmaxp.LT.2) nmaxp=2
176 
177  ! set up local PL, QL
178  allocate(pl_loc(0:nmaxp+1))
179  allocate(ql_loc(0:nmaxp+1))
180 
181  IF (mode_loc.EQ.0) THEN
182  gamma=gammah(m,over)*ar*pisq
183  IF (abs(gamma).LT.tiny) THEN
184  ierr =1
185  err_msg = 'M IS TOO LARGE FOR MODE=0, BETTER TRY MODE=1'
186  chckerr(err_msg)
187  END IF
188  ELSE
189  gamma=ar
190  END IF
191  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192  ! WE EVALUATE THE CONTINUED FRACTION USING !
193  ! LENTZ-THOMPSON !
194  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
195  ierr = frac(z,m,0,eps,tinysq,fc)
196  chckerr('')
197  qdc1=qz*qz-1._dp
198  qargu=qz/dsqrt(qdc1)
199  !DFAC1=DPPI*GAMMA/PI ! Is never used?
200  !DFAC2=GAMMA/DPPI ! Is never used?
201  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
202  ! WE EVALUATE Q_{-1/2},Q^{1}_{-1/2} !
203  ! USING SLATEC ROUTINES FOR ELLIPTIC FUNCTIONS !
204  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
205  argu1=dsqrt(2._dp/(z+1._dp))
206  qlmm(0)=argu1*ellip1(argu1)
207  qlmm(1)=-1._dp/dsqrt(2._dp*(qz-1._dp))*ellip2(argu1)
208  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
209  ! WE APPLY FORWARD RECURRENCE IN M FOR Q'S !
210  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
211  mp=1
212  IF (mode_loc.EQ.0) THEN
213  1 IF ((mp.LE.m).AND.(abs(qlmm(mp)).LT.over)) THEN
214  qlmm(mp+1)=-2._dp*mp*qargu*qlmm(mp)-&
215  &(mp-0.5_dp)*(mp-0.5_dp)*qlmm(mp-1)
216  mp=mp+1
217  GOTO 1
218  ENDIF
219  IF ((mp-1).NE.m) THEN
220  ierr = 1
221  err_msg = 'M IS TOO LARGE FOR MODE=0, BETTER TRY MODE=1'
222  chckerr(err_msg)
223  END IF
224  ELSE
225  qlmm(0)=qlmm(0)/pisq
226  qlmm(1)=qlmm(1)*2._dp/pisq
227  2 IF ((mp.LE.m).AND.(abs(qlmm(mp)).LT.over)) THEN
228  d1=mp+0.5_dp
229  qlmm(mp+1)=-2._dp*mp*qargu*qlmm(mp)/d1&
230  &-(mp-0.5_dp)*qlmm(mp-1)/d1
231  mp=mp+1
232  GOTO 2
233  ENDIF
234  IF ((mp-1).NE.m) THEN
235  ierr = 1
236  err_msg = 'M IS TOO LARGE FOR MODE=1,2'
237  chckerr(err_msg)
238  END IF
239  END IF
240  !NMMAX=M ! Is never used?
241  dfacqs=-(gamma/qlmm(m+1))*gamma/pi
242  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
243  ! EVALUATION OF PL(0) !
244  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
245  IF (ical.EQ.1) THEN
246  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
247  ! WE CALCULATE THE CF FOR P'S !
248  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
249  ierr = fracps(qargu,m+1,0,eps,tinysq,fcp)
250  chckerr('')
251  dd=m+0.5_dp
252  IF (mode_loc.NE.0) THEN
253  fcp=fcp/dd
254  dfacqs=dfacqs/dd
255  END IF
256  pl_loc(0)=dfacqs/dsqrt(qdc1)/(1._dp-fcp*qlmm(m)/qlmm(m+1))
257  ELSE
258  CALL expan(z,mode_loc,ipre_loc,over,qargu,m,pl0)
259  pl_loc(0)=pl0
260  END IF
261  qm0=qlmm(m)
262  dfac3=(gamma/qm0)*gamma/pi/(0.5_dp-m)
263  pl_loc(1)=pl_loc(0)*fc+dfac3
264  np=1
265  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
266  ! WE USE THE RECURRENCE RELATIONS FOR P'S !
267  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
268  3 IF ((np.LE.nmaxp).AND.(abs((np-m+0.5_dp)*pl_loc(np)).LT.over)) THEN
269  pl_loc(np+1)=(2._dp*np*z*pl_loc(np)-(np+m-0.5_dp)*pl_loc(np-1))/&
270  &(np-m+0.5_dp)
271  np=np+1
272  GOTO 3
273  ENDIF
274  nmaxp=np-1
275  newn=nmaxp+1 ! incremented by 1 to get more logical output
276  dfac4=(factco(nmaxp,pl_loc(nmaxp+1),m)*gamma)*gamma/pi/(nmaxp+m+0.5_dp)
277  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
278  ! WE EVALUATE THE CONTINUED FRACTION USING LENTZ-THOMPSON !
279  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280  ierr = frac(z,m,nmaxp,eps,tinysq,fc)
281  chckerr('')
282  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
283  ! EVALUATION OF PL(NMAX+1) AND PL(NMAX) USING !
284  ! THE WRONSKIAN W{PL(NMAX),QL(NMAX)}, !
285  ! THE KNOWN VALUES OF PL(NMAX+1) AND PL(NMAX) !
286  ! THE VALUE OF H = QL(NMAX+1)/QL(NMAX) !
287  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288  ql_loc(nmaxp)=dfac4/(1._dp-fc*pl_loc(nmaxp)/pl_loc(nmaxp+1))
289  ql_loc(nmaxp+1)=ql_loc(nmaxp)*fc
290  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291  ! WE USE THE BACKWARD RECURRENCE RELATION !
292  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293  DO 4 i=1,nmaxp
294  np=nmaxp+1-i
295  n=np-1
296  ql_loc(n)=((np+np)*z*ql_loc(np)-(np-m+0.5_dp)*ql_loc(np+1))/&
297  &(np+m-0.5_dp)
298  4 CONTINUE
299 
300  ! copy local PL and QL back
301  pl = pl_loc(0:nmax)
302  ql = ql_loc(0:nmax)
303  end function dtorh1
304 
305  ! Adapted to 100k iterations in stead of 10k.
306  integer function frac(Z,M,NMAX,EPS,TINYSQ,FC) result(ierr)
307  character(*), parameter :: rout_name = 'FRAC'
308 
309  INTEGER m,nmax,n,mm
310  REAL(dp) z,eps,tinysq,fc,dz2,dn0,dn1,dn2,dn3,dn4,b,a,c0,d0,delta
311 
312  ! initialize ierr
313  ierr = 0
314 
315  n=nmax
316  mm=0
317  dz2=z+z
318  dn0=n+m
319  dn1=n+1._dp
320  dn2=dn0+0.5_dp
321  dn3=dn0-0.5_dp
322  dn4=1._dp-m-m
323  b=2._dp*dn1*z/dn2
324  a=1._dp
325  fc=tinysq
326  c0=fc
327  d0=0._dp
328  81 d0=b+a*d0
329  IF (dabs(d0).LT.tinysq) d0=tinysq
330  c0=b+a/c0
331  IF (dabs(c0).LT.tinysq) c0=tinysq
332  d0=1._dp/d0
333  delta=c0*d0
334  fc=fc*delta
335  mm=mm+1
336  a=-(1.d0+dn4/(dn3+mm))
337  b=dz2*(dn1+mm)/(dn2+mm)
338  IF (mm.LT.1000000) THEN
339  IF (dabs(delta-1.d0).GT.eps) GOTO 81
340  END IF
341  IF (mm.EQ.1000000) then
342  ierr =1
343  err_msg = 'CF CONVERGENCE FAILS FOR Z = '//trim(r2str(z))
344  chckerr(err_msg)
345  END IF
346  end function frac
347 
348  integer function fracps(QZ,M,N,EPS,TINYSQ,FC) result(ierr)
349  character(*), parameter :: rout_name = 'FRACPS'
350 
351  INTEGER m,n,mm
352  REAL(dp) qz,eps,tinysq,fc,dn2,dz2,dm,b,a,c0,d0,delta
353 
354  ! initialize ierr
355  ierr = 0
356 
357  mm=0
358  dn2=n*n
359  dz2=qz+qz
360  dm=m-0.5_dp
361  b=dz2*m
362  a=dn2-dm*dm
363  fc=tinysq
364  c0=fc
365  d0=0._dp
366  82 d0=b+a*d0
367  IF (dabs(d0).LT.tinysq) d0=tinysq
368  c0=b+a/c0
369  IF (dabs(c0).LT.tinysq) c0=tinysq
370  d0=1._dp/d0
371  delta=c0*d0
372  fc=fc*delta
373  mm=mm+1
374  a=dn2-(mm+dm)*(mm+dm)
375  b=dz2*(mm+m)
376  IF (mm.LT.1000000) THEN
377  IF (abs(delta-1._dp).GT.eps) GOTO 82
378  END IF
379  IF (mm.EQ.1000000) then
380  ierr =1
381  err_msg = 'CF CONVERGENCE FAILS FOR QZ = '//trim(r2str(qz))
382  chckerr(err_msg)
383  END IF
384  end function fracps
385 
386  REAL(dp) function ellip2(ak)
387  REAL(dp) ak
388  REAL(dp) q
389 
390  q=(1._dp-ak)*(1._dp+ak)
391  ellip2=drf(q)-(ak)**2*drd(q)/3.d0
392  END FUNCTION ellip2
393 
394  REAL(dp) function ellip1(ak)
395  REAL(dp) ak
396 
397  ellip1=drf((1._dp-ak)*(1._dp+ak))
398  END FUNCTION ellip1
399 
400 
401  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
402  ! WE USE SLATEC ROUTINES IN THE EVALUATION OF ELLIPTIC INTEGRALS !
403  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
404  REAL(dp) function drf (y)
405  REAL(dp) epslon, errtol
406  REAL(dp) c1, c2, c3, e2, e3, lamda
407  REAL(dp) mu, s, xn, xndev
408  REAL(dp) xnroot, y, yn, yndev, ynroot, zn, zndev, znroot
409  LOGICAL first
410  SAVE errtol,c1,c2,c3,first
411  DATA first /.true./
412 
413  !***FIRST EXECUTABLE STATEMENT DRF
414  IF (first) THEN
415  errtol = 1.e-8_dp
416  c1 = 1.0_dp/24.0_dp
417  c2 = 3.0_dp/44.0_dp
418  c3 = 1.0_dp/14.0_dp
419  ENDIF
420  first = .false.
421  drf = 0.0_dp
422  xn = 0._dp
423  yn = y
424  zn = 1._dp
425 
426  30 mu = (xn+yn+zn)/3.0_dp
427  xndev = 2.0_dp - (mu+xn)/mu
428  yndev = 2.0_dp - (mu+yn)/mu
429  zndev = 2.0_dp - (mu+zn)/mu
430  epslon = max(abs(xndev),abs(yndev),abs(zndev))
431  IF (epslon.LT.errtol) GO TO 40
432  xnroot = sqrt(xn)
433  ynroot = sqrt(yn)
434  znroot = sqrt(zn)
435  lamda = xnroot*(ynroot+znroot) + ynroot*znroot
436  xn = (xn+lamda)*0.250_dp
437  yn = (yn+lamda)*0.250_dp
438  zn = (zn+lamda)*0.250_dp
439  GO TO 30
440 
441  40 e2 = xndev*yndev - zndev*zndev
442  e3 = xndev*yndev*zndev
443  s = 1.0_dp + (c1*e2-0.10_dp-c2*e3)*e2 + c3*e3
444  drf = s/sqrt(mu)
445  END FUNCTION drf
446 
447 
448  REAL(dp) function drd (y)
449  REAL(dp) epslon, errtol
450  REAL(dp) c1, c2, c3, c4, ea, eb, ec, ed, ef, lamda
451  REAL(dp) mu, power4, sigma, s1, s2, xn, xndev
452  REAL(dp) xnroot, y, yn, yndev, ynroot, zn, zndev, znroot
453  LOGICAL first
454  SAVE errtol, c1, c2, c3, c4, first
455  DATA first /.true./
456 
457  !***FIRST EXECUTABLE STATEMENT DRD
458  IF (first) THEN
459  errtol = 1.e-8_dp
460  c1 = 3.0_dp/14.0_dp
461  c2 = 1.0_dp/6.0_dp
462  c3 = 9.0_dp/22.0_dp
463  c4 = 3.0_dp/26.0_dp
464  ENDIF
465  first = .false.
466 
467  ! CALL ERROR HANDLER IF NECESSARY.
468 
469  drd = 0.0_dp
470 
471  xn = 0._dp
472  yn = y
473  zn = 1._dp
474  sigma = 0.0_dp
475  power4 = 1.0_dp
476 
477  30 mu = (xn+yn+3.0_dp*zn)*0.20_dp
478  xndev = (mu-xn)/mu
479  yndev = (mu-yn)/mu
480  zndev = (mu-zn)/mu
481  epslon = max(abs(xndev), abs(yndev), abs(zndev))
482  IF (epslon.LT.errtol) GO TO 40
483  xnroot = sqrt(xn)
484  ynroot = sqrt(yn)
485  znroot = sqrt(zn)
486  lamda = xnroot*(ynroot+znroot) + ynroot*znroot
487  sigma = sigma + power4/(znroot*(zn+lamda))
488  power4 = power4*0.250_dp
489  xn = (xn+lamda)*0.250_dp
490  yn = (yn+lamda)*0.250_dp
491  zn = (zn+lamda)*0.250_dp
492  GO TO 30
493 
494  40 ea = xndev*yndev
495  eb = zndev*zndev
496  ec = ea - eb
497  ed = ea - 6.0_dp*eb
498  ef = ed + ec + ec
499  s1 = ed*(-c1+0.250_dp*c3*ed-1.50_dp*c4*zndev*ef)
500  s2 = zndev*(c2*ef+zndev*(-c3*ec+zndev*c4*ea))
501  drd = 3.0_dp*sigma + power4*(1.0_dp+s1+s2)/(mu*sqrt(mu))
502  END FUNCTION drd
503 
504  REAL(dp) function gammah(n,over)
505  INTEGER n,i,j
506  REAL(dp) over
507 
508  i=n
509  j=2*i-1
510  gammah=1._dp
511  85 IF ((j.GE.1).AND.(gammah.LT.over)) THEN
512  gammah=gammah*j/2._dp
513  i=i-1
514  j=2*i-1
515  GOTO 85
516  ENDIF
517  IF (j.GT.1) THEN
518  gammah=0._dp
519  END IF
520  END FUNCTION gammah
521 
522  REAL(dp) function factco(n,pl,m)
523  INTEGER n,m,j
524  REAL(dp) pl,x1,x2
525 
526  factco=1._dp/pl
527  x1=m+0.5_dp
528  x2=-m+0.5_dp
529  j=n
530  86 IF (j.GE.0) THEN
531  factco=factco*(j+x1)/(j+x2)
532  j=j-1
533  GOTO 86
534  ENDIF
535  END FUNCTION factco
536 
537  SUBROUTINE expan(Z,MODE,IPRE,OVER,QARGU,M,PL)
538  INTEGER mode,ipre,m,k,i
539  REAL(dp) z,over,qargu,pl,pisq,db,fl,cc,dz,gamma,ar,&
540  &DFAC,DF1,A0,Z2I,DA1,DA2,DELTA,SUM,DCC,DCCP,DKK
541  REAL(dp) preci(2)
542 
543  preci(1)=1.e-13_dp
544  preci(2)=1.e-9_dp
545  pisq=dsqrt(pi)
546  db=2._dp*dlog(2._dp)
547  fl=m/2.
548  cc=abs(float(int(fl))-fl)
549  IF (cc.LT.0.4_dp) THEN
550  ar=1.
551  ELSE
552  ar=-1.
553  END IF
554  dz=1._dp
555  DO 100 i=1,m
556  dz=dz/qargu
557  100 CONTINUE
558  IF (mode.EQ.0) THEN
559  gamma=gammah(m,over)*ar*pisq
560  ELSE
561  gamma=ar
562  END IF
563  dfac=2._dp/pi*dz*gamma/pisq
564  df1=dlog(2._dp*z)
565  a0=1._dp/dsqrt(2._dp*z)
566  z2i=1._dp/(z*z)
567  delta=1._dp
568  sum=0._dp
569  k=0
570  da2=factor(m,0)
571  da1=db+psi(m,0)
572  87 delta=(df1+da1)*da2*a0
573  sum=sum+delta
574  dcc=0.5_dp+m+2._dp*k
575  dccp=dcc+1._dp
576  dkk=k+1._dp
577  da2=da2*dccp*dcc/(dkk*dkk)
578  da1=da1+1._dp/dkk-1._dp/dccp-1._dp/dcc
579  k=k+1
580  a0=a0*0.25_dp*z2i
581  IF (abs(delta/sum).GT.preci(ipre)) GOTO 87
582  pl=sum*dfac
583  END SUBROUTINE expan
584 
585  REAL(dp) function factor(m,k)
586  INTEGER m,k,n,i,j
587  REAL(dp) x1
588 
589  factor=1._dp
590  IF (k.GE.1) THEN
591  x1=m+0.5_dp
592  n=2*k-1
593  i=k
594  j=n
595  88 IF (j.GE.0) THEN
596  factor=factor*(j+x1)/i/i
597  j=j-1
598  i=i-1
599  IF (i.EQ.0) i=1
600  GOTO 88
601  ENDIF
602  END IF
603  END FUNCTION factor
604 
605  REAL(dp) function psi(m,k)
606  INTEGER m,k,n,j,i
607  REAL(dp) factr1,factr2
608  psi=0._dp
609  factr1=0._dp
610  factr2=0._dp
611  n=2*k+m
612  j=1
613  IF (k.GE.1) THEN
614  89 IF (j.LE.k) THEN
615  factr1=factr1+1._dp/j
616  j=j+1
617  GOTO 89
618  ENDIF
619  END IF
620  i=1
621  90 IF (i.LE.n) THEN
622  factr2=factr2+1._dp/(2._dp*i-1._dp)
623  i=i+1
624  GOTO 90
625  ENDIF
626  IF (n.EQ.0) factr2=0._dp
627  psi=factr1-2._dp*factr2
628  END FUNCTION psi
629 end module
dtorh
Calculation of toroidal functions and .
Definition: dtorh.f90:8
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
dtorh::dtorh1
integer function, public dtorh1(z, m, nmax, pl, ql, newn, mode, ipre)
Calculates toroidal harmonics of a fixed order m and degrees up to nmax.
Definition: dtorh.f90:75
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
num_vars::pi
real(dp), parameter, public pi
Definition: num_vars.f90:83
this documentation was created with Doxygen 1.8.17 by Toon Weyens