Newer
Older
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
139
140
141
142
143
144
145
!! This source has been adapted from the GENE code https://genecode.org/ !!
!>Provides a circular concentric geometry
MODULE circular
USE prec_const
USE basic
USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0
! use coordinates,only: gcoor, get_dzprimedz
USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven
! use discretization
USE lagrange_interpolation
implicit none
public:: get_circ
private
! SUBROUTINES DEFINITIONS
CONTAINS
SUBROUTINE get_circ(trpeps,major_R,major_Z,q0,shat,&
C_y,C_xy,Cyq0_x0,gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,dBdx_,dBdy_,dBdz_,&
Bfield_,jacobian_,R_hat_,Z_hat_,dxdR_,dxdZ_)
!!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
real(xp), INTENT(INOUT) :: major_R ! major radius
real(xp), INTENT(INOUT) :: major_Z ! major Z
real(xp), INTENT(INOUT) :: q0 ! safetyfactor
real(xp), INTENT(INOUT) :: shat ! safetyfactor
real(xp), INTENT(INOUT) :: C_y, C_xy, Cyq0_x0
real(xp), dimension(local_nz+ngz,nzgrid), INTENT(OUT) :: &
gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,&
dBdx_,dBdy_,dBdz_,Bfield_,jacobian_,&
R_hat_,Z_hat_,dxdR_,dxdZ_
INTEGER :: iz,eo
! No parameter in gyacomo yet
real(xp) :: SIGN_Ip=1._xp ! current sign (only normal current)
!!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
REAL(xp) :: g11,g12,g22,g33 ! normalized metric in (PSI,CHI,PHI)
REAL(xp) :: B,dBdr_c,dBdchi ! J_PCP : J_PSI,CHI,P
! intermediate variables, t stands for theta the poloidal angle.
REAL(xp) :: qbar,dpsidr,dchidr,dchidt,dBdr,dBdt,cost,sint
REAL(xp) :: fac2, dCydr_Cy
REAL(xp) :: dxdr, dqdr, dqdx
REAL(xp) :: dqbardr, minor_r, r0, z
minor_r = major_R
r0 = major_R*trpeps
qbar = abs(q0)*sqrt(1-trpeps**2)
C_y = abs(r0/q0)*SIGN_Ip
dxdr = 1._xp ! x=r
dpsidr = r0/qbar ! 1/(Bref*Lref)*dpsi/dr
! dqdx = shat*q0/minor_r/dxdr ! AH: added
dqdx = shat/C_y ! AH: added
dqdr = dqdx*dxdr/minor_r
dqbardr = sqrt(1._xp-trpeps**2)*(abs(dqdr)-abs(q0)*trpeps/major_R/(1-trpeps**2))
dCydr_Cy = 0._xp
C_xy = dpsidr/(dxdr*abs(C_y))
Cyq0_x0 = C_y*q0/r0
fac2 = (dCydr_Cy+dqdr/q0)
DO eo = 1,nzgrid
do iz=1,local_nz+ngz
z = zarray(iz,eo)
! z dep. intermediate variables
cost = (cos(z) - trpeps)/(1._xp-trpeps*cos(z))
sint = sqrt(1._xp-trpeps**2)*sin(SIGN_Ip*z)/(1._xp-trpeps*cos(z))
! derivatives with respect to r,theta
dchidr = -1._xp/major_R*sin(z)/(1._xp-trpeps**2)
dchidt = SIGN_Ip * sqrt(1-trpeps**2)/(1+trpeps*cost)
B = 1._xp/(1._xp+trpeps*cost)*sqrt(1.+(trpeps/qbar)**2)
dBdr = B/major_R*(-cost/(1._xp+trpeps*cost) + &
+ trpeps/(trpeps**2+qbar**2)*(1._xp - major_R*trpeps*dqbardr/qbar))
dBdt = trpeps*sint*B / (1._xp+trpeps*cost)
! metric in (r,CHI,PHI)
g11 = 1._xp
g22 = dchidr**2+1._xp/(major_R*trpeps)**2*dchidt**2
g12 = dchidr
g33 = 1._xp/major_R**2*1._xp/(1._xp+trpeps*cost)**2
! magnetic field derivatives in
dBdchi = dBdt/dchidt
dBdr_c = 1._xp/g11*(dBdr-g12*dBdt/dchidt)
! metric in (x,y,z)
gxx_(iz,eo) = (dxdr)**2*g11
gyy_(iz,eo) = (C_y*q0)**2 * ((fac2*z)**2*g11 +2._xp*fac2*z*g12+g22) + C_y**2*g33
gxy_(iz,eo) = dxdr*C_y*SIGN_Ip*q0*(fac2*z*g11+g12)
gxz_(iz,eo) = dxdr*g12
gyz_(iz,eo) = C_y*q0*SIGN_Ip*(fac2*z*g12+g22)
gzz_(iz,eo) = g22
! jacobian
jacobian_(iz,eo)=C_xy*abs(q0)*major_R*(1+trpeps*cost)**2
! Bfield
Bfield_(iz,eo) = B
dBdx_(iz,eo) = dBdr_c/dxdr
dBdy_(iz,eo) = 0._xp
dBdz_(iz,eo) = dBdchi
!derivatives with respect to cylindrical coordinates
dxdR_(iz,eo) = cost
dxdZ_(iz,eo) = sint
R_hat_(iz,eo) = major_R*(1.+trpeps*cost)
Z_hat_(iz,eo) = major_Z+major_R*trpeps*sint
!!!!!!! TEST salpha
! ! metric
! gxx_(iz,eo) = 1._xp
! gxy_(iz,eo) = shat*z - amhd*SIN(z)
! gxz_(iz,eo) = 0._xp
! gyy_(iz,eo) = 1._xp + (shat*z - amhd*SIN(z))**2
! gyz_(iz,eo) = 1._xp/trpeps
! gzz_(iz,eo) = 1._xp/trpeps**2
! dxdR_(iz,eo)= COS(z)
! dxdZ_(iz,eo)= SIN(z)
! ! ! Poloidal plane coordinates
! R_hat_(iz,eo) = 1._xp + trpeps*COS(z)
! Z_hat_(iz,eo) = trpeps*SIN(z)
! ! Relative strengh of modulus of B
! Bfield_(iz,eo) = 1._xp/(1._xp + trpeps*COS(z))
! ! Jacobian
! Jacobian_(iz,eo) = q0/Bfield_(iz,eo)
! ! Derivative of the magnetic field strenght
! dBdx_(iz,eo) = -COS(z)*Bfield_(iz,eo)**2 ! LB = 1
! dBdy_(iz,eo) = 0._xp
! dBdz_(iz,eo) = trpeps*SIN(z)*Bfield_(iz,eo)**2
! ! Curvature factor
! C_xy = 1._xp
!!!!! END TEST
end do
END DO
END SUBROUTINE get_circ
!-----------------------------------------------------------
END MODULE circular