forked from PrincetonUniversity/SPEC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstzxyz.f90
138 lines (103 loc) · 5.55 KB
/
stzxyz.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
!> \file
!> \brief Calculates coordinates, \f${\bf x}(s,\theta,\zeta) \equiv R \, {\bf e}_R + Z \, {\bf e}_Z\f$, and metrics, at given \f$(s,\theta,\zeta)\f$.
!> \brief Calculates coordinates, \f${\bf x}(s,\theta,\zeta) \equiv R \, {\bf e}_R + Z \, {\bf e}_Z\f$, and metrics, at given \f$(s,\theta,\zeta)\f$.
!> \ingroup grp_diagnostics
!>
!> <ul>
!> <li> This routine is a "copy" of co01aa(),
!> which calculates the coordinate information on a regular, discrete grid in \f$\theta\f$ and \f$\zeta\f$ at given \f$s\f$
!> whereas stzxyz() calculates the coordinate information at a single point \f$(s,\theta,\zeta)\f$. </li>
!> <li> \todo Please see co01aa() for documentation.
!>
!> </li>
!> </ul>
!>
!> @param[in] lvol
!> @param[in] stz
!> @param[out] RpZ
subroutine stzxyz( lvol , stz , RpZ )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
use constants, only : zero, one, half
use numerical, only : vsmall
use fileunits, only : ounit
use inputlist, only : Wstzxyz, Igeometry, Nvol, Ntor
use cputiming, only : Tstzxyz
use allglobal, only : myid, cpus, mn, im, in, halfmm, MPI_COMM_SPEC, &
iRbc, iZbs, iRbs, iZbc, &
Lcoordinatesingularity, &
NOTstellsym, &
Mvol
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
LOCALS
INTEGER, intent(in) :: lvol
REAL, intent(in) :: stz(1:3)
REAL, intent(out) :: RpZ(1:3)
INTEGER :: ii, mi, ni
REAL :: Remn, Zomn, Romn, Zemn, RR, phi, ZZ, arg, carg, sarg, lss, alss, blss, sbar, sbarhim, fj
BEGIN(stzxyz)
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
#ifdef DEBUG
FATAL(stzxyz, lvol.lt.1 .or. lvol.gt.Mvol, invalid interface label )
FATAL(stzxyz, abs(stz(1)).gt.one, invalid radial coordinate )
#endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
RpZ(1:3) = zero ; RR = zero ; phi = stz(3) ; ZZ = zero ! initialize intent(out), summations ; 17 Dec 15;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
lss = stz(1) ; alss = half - lss * half ; blss = half + lss * half ! shorthand; 17 Dec 15;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
do ii = 1, mn ; mi = im(ii) ; ni = in(ii) ; arg = mi * stz(2) - ni * phi ; carg = cos(arg) ; sarg = sin(arg) ! Fourier sum; 17 Dec 15;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Remn = zero ; Zomn = zero
! if( NOTstellsym ) then
Romn = zero ; Zemn = zero
! endif
if( Lcoordinatesingularity ) then
sbar = ( lss + one ) * half
if( mi.eq.0 ) then
if( Igeometry.eq.2 ) then ; fj = sbar
else ; fj = sbar**2
endif
else ; fj = sbar**im(ii)
endif
Remn = iRbc(ii,0) + ( iRbc(ii,1) - iRbc(ii,0) ) * fj
if( NOTstellsym ) then
Romn = iRbs(ii,0) + ( iRbs(ii,1) - iRbs(ii,0) ) * fj
endif
if( Igeometry.eq.3 ) then ! recall that for cylindrical geometry there is no need for Z; 20 Apr 13;
Zomn = iZbs(ii,0) + ( iZbs(ii,1) - iZbs(ii,0) ) * fj
if( NOTstellsym ) then
Zemn = iZbc(ii,0) + ( iZbc(ii,1) - iZbc(ii,0) ) * fj
endif
endif
else ! matches if( Lcoordinatesingularity) ; 20 Apr 13;
Remn = alss * iRbc(ii,lvol-1) + blss * iRbc(ii,lvol)
if( NOTstellsym ) then
Romn = alss * iRbs(ii,lvol-1) + blss * iRbs(ii,lvol)
else
Romn = zero
endif ! end of if( NOTstellsym ) ; 22 Apr 13;
if( Igeometry.eq.3 ) then
Zomn = alss * iZbs(ii,lvol-1) + blss * iZbs(ii,lvol)
if( NOTstellsym ) then
Zemn = alss * iZbc(ii,lvol-1) + blss * iZbc(ii,lvol)
else
Zemn = zero
endif ! end of if( NOTstellsym ) ; 22 Apr 13;
else
Zomn = zero
Zemn = zero
endif ! end of if( Igeometry.eq.3 ) ; 22 Apr 13;
endif ! end of if( Lcoordinatesingularity ) ; 20 Feb 13;
RR = RR + Remn * carg + Romn * sarg
if( Igeometry.eq.3 ) then
ZZ = ZZ + Zemn * carg + Zomn * sarg
endif
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
enddo ! end of do ii = 1, mn ; Fourier summation loop;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
RpZ(1:3) = (/ RR, phi, ZZ /) ! return coordinate functions;
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
RETURN(stzxyz)
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
end subroutine stzxyz
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!