From c3e24483baac8eb5aa0bad90de620e3f03d97ab4 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 14 Oct 2022 16:59:14 +0200
Subject: [PATCH] Miller geometry is implemented, not benchmarked yet
---
src/geometry_mod.F90 | 45 +-
src/{lag_interp.F90 => lag_interp_mod.F90} | 60 +-
src/miller_geometry.F90 | 589 -------------------
src/miller_mod.F90 | 645 +++++++++++++++++++++
4 files changed, 710 insertions(+), 629 deletions(-)
rename src/{lag_interp.F90 => lag_interp_mod.F90} (78%)
delete mode 100644 src/miller_geometry.F90
create mode 100644 src/miller_mod.F90
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 665f18ba..9840cc0b 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -9,17 +9,34 @@ module geometry
use fields
use basic
use calculus, ONLY: simpson_rule_z
-
+ use miller, ONLY: set_miller_parameters, get_miller
implicit none
PRIVATE
- ! Geometry parameters
+ ! Geometry input parameters
CHARACTER(len=16), &
PUBLIC, PROTECTED :: geom
- REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
- REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
- REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
- LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
- ! Geometrical operators
+ REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
+ REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
+ REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
+ REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
+ ! parameters for Miller geometry
+ REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation
+ REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr
+ REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity
+ REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
+ REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness
+ REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr
+
+ ! GENE unused additional parameters for miller_mod
+ REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0 ! meant to redistribute the points in z
+ REAL(dp), PUBLIC, PROTECTED :: major_R = 1 ! major radius
+ REAL(dp), PUBLIC, PROTECTED :: major_Z = 0 ! vertical elevation
+ REAL(dp), PUBLIC, PROTECTED :: dpdx_pm_geom = 0 ! amplitude mag. eq. pressure grad.
+ REAL(dp), PUBLIC, PROTECTED :: C_y = 0 ! defines y coordinate : Cy (q theta - phi)
+ REAL(dp), PUBLIC, PROTECTED :: C_xy = 0 ! defines x coordinate : B = Cxy Vx x Vy
+
+ ! Geometrical auxiliary variables
+ LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
! Curvature
REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
@@ -38,6 +55,7 @@ implicit none
REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition
INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R
+
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
eval_magnetic_geometry, set_ikx_zBC_map
@@ -47,7 +65,8 @@ CONTAINS
SUBROUTINE geometry_readinputs
! Read the input parameters
IMPLICIT NONE
- NAMELIST /GEOMETRY/ geom, q0, shear, eps
+ NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
+ kappa, s_kappa,delta, s_delta, zeta, s_zeta ! For miller
READ(lu_in,geometry)
IF(shear .NE. 0._dp) SHEARED = .true.
@@ -78,8 +97,14 @@ CONTAINS
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
SHEARED = .FALSE.
+ CASE('miller')
+ IF( my_id .eq. 0 ) WRITE(*,*) 'Miller geometry'
+ call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
+ call get_miller(eps,major_R,major_Z,q0,shear,alpha_MHD,edge_opt,&
+ C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
+ gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ)
CASE DEFAULT
- ERROR STOP 'Error stop: geometry not recognized!!'
+ STOP 'geometry not recognized!!'
END SELECT
ENDIF
!
@@ -116,7 +141,7 @@ CONTAINS
SUBROUTINE eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
- REAL(dp) :: z, kx, ky, alpha_MHD, Gx, Gy
+ REAL(dp) :: z, kx, ky, Gx, Gy
alpha_MHD = 0._dp
parity: DO eo = 0,1
diff --git a/src/lag_interp.F90 b/src/lag_interp_mod.F90
similarity index 78%
rename from src/lag_interp.F90
rename to src/lag_interp_mod.F90
index db2afc48..8415df5d 100644
--- a/src/lag_interp.F90
+++ b/src/lag_interp_mod.F90
@@ -1,8 +1,8 @@
-#include "redef.h"
!! This source is taken from GENE https://genecode.org/ !!
!>lagrange_interpolation contains subroutines to perform
!!a mid-point lagrange interpolation of order 3
MODULE lagrange_interpolation
+ USE prec_const
IMPLICIT NONE
PUBLIC:: lag3interp, lag3deriv, lag3interp_2d
PUBLIC:: lag3interp_complex
@@ -21,11 +21,11 @@ CONTAINS
!> Third order lagrange interpolation
SUBROUTINE lag3interp_scalar(y_in,x_in,n_in,y_out,x_out)
INTEGER, INTENT(IN) :: n_in
- REAL, DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL, INTENT(IN) :: x_out
- REAL, INTENT(OUT) :: y_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), INTENT(IN) :: x_out
+ REAL(dp), INTENT(OUT) :: y_out
- REAL, DIMENSION(1) :: xout_wrap, yout_wrap
+ REAL(dp), DIMENSION(1) :: xout_wrap, yout_wrap
xout_wrap = x_out
call lag3interp_array(y_in,x_in,n_in,yout_wrap,xout_wrap,1)
@@ -36,11 +36,11 @@ CONTAINS
!> Third order lagrange interpolation
subroutine lag3interp_array(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
- REAL, DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL, DIMENSION(n_out), INTENT(IN) :: x_out
- REAL, DIMENSION(n_out), INTENT(OUT) :: y_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(dp), DIMENSION(n_out), INTENT(OUT) :: y_out
- REAL :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
@@ -91,11 +91,11 @@ CONTAINS
SUBROUTINE lag3interp_complex(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
COMPLEX, DIMENSION(n_in), INTENT(IN) :: y_in
- REAL, DIMENSION(n_in), INTENT(IN) :: x_in
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: x_in
COMPLEX, DIMENSION(n_out), INTENT(OUT) :: y_out
- REAL, DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
- REAL :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
@@ -143,22 +143,22 @@ CONTAINS
!>2D interpolation
- !\TODO check whether a "real" 2D interpolation would
+ !\TODO check whether a "REAL(dp)" 2D interpolation would
!! be more appropriate
SUBROUTINE lag3interp_2d(y_in,x1_in,n1_in,x2_in,n2_in,&
&y_out,x1_out,n1_out,x2_out,n2_out)
INTEGER, INTENT(IN) :: n1_in,n2_in,n1_out,n2_out
- REAL, DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
- REAL, DIMENSION(n1_in) :: x1_in
- REAL, DIMENSION(n2_in) :: x2_in
- REAL, DIMENSION(n1_out), INTENT(IN) :: x1_out
- REAL, DIMENSION(n2_out), INTENT(IN) :: x2_out
- REAL, DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
+ REAL(dp), DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
+ REAL(dp), DIMENSION(n1_in) :: x1_in
+ REAL(dp), DIMENSION(n2_in) :: x2_in
+ REAL(dp), DIMENSION(n1_out), INTENT(IN) :: x1_out
+ REAL(dp), DIMENSION(n2_out), INTENT(IN) :: x2_out
+ REAL(dp), DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
!local variables
- REAL, DIMENSION(n2_in) :: y2_in_tmp
- REAL, DIMENSION(n2_out) :: y2_out_tmp
- REAL, DIMENSION(n1_in,n2_out) :: y_tmp
+ REAL(dp), DIMENSION(n2_in) :: y2_in_tmp
+ REAL(dp), DIMENSION(n2_out) :: y2_out_tmp
+ REAL(dp), DIMENSION(n1_in,n2_out) :: y_tmp
INTEGER :: i
DO i=1,n1_in
@@ -181,11 +181,11 @@ CONTAINS
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_in
- REAL, DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL, INTENT(IN) :: x_out
- REAL, INTENT(OUT) :: dydx_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), INTENT(IN) :: x_out
+ REAL(dp), INTENT(OUT) :: dydx_out
- REAL, DIMENSION(1) :: xout_wrap, dydxout_wrap
+ REAL(dp), DIMENSION(1) :: xout_wrap, dydxout_wrap
xout_wrap = x_out
call lag3deriv_array(y_in,x_in,n_in,dydxout_wrap,xout_wrap,1)
@@ -198,11 +198,11 @@ CONTAINS
!>Returns Derivative based on a 3rd order lagrange interpolation
subroutine lag3deriv_array(y_in,x_in,n_in,dydx_out,x_out,n_out)
INTEGER :: n_in,n_out
- REAL, DIMENSION(n_in), INTENT(IN) :: y_in,x_in
- REAL, DIMENSION(n_out), INTENT(IN) :: x_out
- REAL, DIMENSION(n_out), INTENT(OUT) :: dydx_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(dp), DIMENSION(n_out), INTENT(OUT) :: dydx_out
- REAL :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
diff --git a/src/miller_geometry.F90 b/src/miller_geometry.F90
deleted file mode 100644
index 15f70539..00000000
--- a/src/miller_geometry.F90
+++ /dev/null
@@ -1,589 +0,0 @@
-#include "redef.h"
-!! This source has been adapted from GENE https://genecode.org/ !!
-!>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998)
-!>and [J. Candy, PPCF 51, 105009 (2009)]
-MODULE miller_mod
- ! use coordinates,only: gcoor, get_dzprimedz
- ! use discretization
- use lagrange_interpolation
- ! use par_geom
- ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, n_pol
- ! use par_other, only: print_ini_msg
-
- implicit none
- public:: get_miller, set_miller_defaults
- public:: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
- public:: thetaShift
- public:: mMode, nMode
- public:: thetak, thetad
- public:: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
- public:: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
- public:: cN_m, sN_m, cNdr_m, sNdr_m
-
- private
-
- REAL :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
-
- INTEGER :: mMode, nMode
- REAL :: thetaShift
- REAL :: thetak, thetad
- REAL :: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
- REAL :: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
-
- INTEGER, PARAMETER :: IND_M=32
- REAL, DIMENSION(0:IND_M-1) :: cN_m, sN_m, cNdr_m, sNdr_m
-
-CONTAINS
-
- !>Set defaults for miller parameters
- subroutine set_miller_defaults
- rho = -1.0
- kappa = 1.0
- s_kappa = 0.0
- delta = 0.0
- s_delta = 0.0
- drR = 0.0
- drZ = 0.0
- zeta = 0.0
- s_zeta = 0.0
-
-
- thetak = 0.0
- thetad = 0.0
-
- aSurf = 0.54
- Delta2 = 1.0
- Delta3 = 1.0
- theta2 = 0.0
- theta3 = 0.0
- Raxis = 1.0
- Zaxis = 0.0
-
- mMode = 2
- nMode = 3
- Deltam = 1.0
- Deltan = 1.0
- s_Deltam = 0.0
- s_Deltan = 0.0
- thetam = 0.0
- thetan = 0.0
-
- cN_m = 0.0
- sN_m = 0.0
- cNdr_m = 0.0
- sNdr_m = 0.0
- end subroutine set_miller_defaults
-
- !>Get Miller metric, magnetic field, jacobian etc.
- subroutine get_miller(geom,edge_opt) !previously: get_miller_a
- type(geomtype),intent(inout):: geom
- real, intent(in):: edge_opt
- integer:: np, np_s, n_pol_ext, n_pol_s
-
- real, dimension(500*(n_pol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi
- real, dimension(500*(n_pol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt
- real, dimension(500*(n_pol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt
- real, dimension(500*(n_pol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
- real, dimension(500*(n_pol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
- real, dimension(500*(n_pol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
- real, dimension(500*(n_pol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
- real, dimension(500*(n_pol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
- real, dimension(500*(n_pol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
- real, dimension(500*(n_pol+2)):: tmp_reverse, theta_reverse, tmp_arr
-
- real, dimension(500*(n_pol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse
- real, dimension(500*(n_pol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s
- real, dimension(500*(n_pol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s
- real, dimension(500*(n_pol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s
- real, dimension(500*(n_pol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s, dtrShape_s
- real, dimension(500*(n_pol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s
- real, dimension(500*(n_pol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s
- real, dimension(500*(n_pol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
- real, dimension(500*(n_pol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
-
- real, dimension(0:nz0-1):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
- real, dimension(0:nz0-1):: R_out, Z_out, dxdR_out, dxdZ_out
- real:: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
- !real :: Lnorm, Psi0 ! currently module-wide defined anyway
- real:: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
- real:: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
- logical:: bMaxShift
- integer:: i, k, iBmax
-
- n_pol_ext = n_pol+2
- n_pol_s = n_pol+1
- np = 500*n_pol_ext
- np_s = 500*n_pol_s
-
- if (rho.lt.0.0) rho = trpeps*major_R
- if (rho.le.0.0) stop 'flux surface radius not defined'
- trpeps = rho/major_R
- gcoor%x0=rho
-
- q0 = sign_Ip_CW * sign_Bt_CW * abs(q0)
-
- R0=major_R
- B0=1.0*sign_Bt_CW
- F=R0*B0
- Z0=major_Z
- pi = acos(-1.0)
- mu_0=4.0*pi
-
- theta=linspace(-pi*n_pol_ext,pi*n_pol_ext-2*pi*n_pol_ext/np,np)
- d_inv=asin(delta)
-
- thetaShift = 0.0
- iBmax = 1
-
- !flux surface parametrization
- thAdj = theta + thetaShift
- select case (magn_geometry)
- case ('miller')
- if (zeta/=0.0 .or. s_zeta/=0.0) then
- R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
- Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
-
- R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
- Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
- + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
-
- R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
- Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
-
- R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
- + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
- Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
- - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
- else
- Rcirc = rho*Cos(thAdj - thetad + thetak)
- Zcirc = rho*Sin(thAdj - thetad + thetak)
- Relong = Rcirc
- Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
- RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
- ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
- Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
- Ztri = ZelongTilt
- RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
- ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
- R = R0 + RtriTilt
- Z = Z0 + ZtriTilt
-
- drRcirc = Cos(thAdj - thetad + thetak)
- drZcirc = Sin(thAdj - thetad + thetak)
- drRelong = drRcirc
- drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
- drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
- drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
- drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
- - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
- drZtri = drZelongTilt
- drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
- drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
- R_rho = drR + drRtriTilt
- Z_rho = drZ + drZtriTilt
-
- dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
- dtZcirc = rho*Cos(thAdj - thetad + thetak)
- dtRelong = dtRcirc
- dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
- dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
- dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
- dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
- dtZtri = dtZelongTilt
- dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
- dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
- R_theta = dtRtriTilt
- Z_theta = dtZtriTilt
-
- dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
- dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
- dtdtRelong = dtdtRcirc
- dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
- dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
- dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
- dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
- + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
- dtdtZtri = dtdtZelongTilt
- dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
- dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
- R_theta_theta = dtdtRtriTilt
- Z_theta_theta = dtdtZtriTilt
- endif
- case default
- write (*,*) "ERROR: invalid analytic geometry specification"
- end select
-
- !dl/dtheta
- dlp=(R_theta**2+Z_theta**2)**0.5
-
- !curvature radius
- Rc=dlp**3*(R_theta*Z_theta_theta-Z_theta*R_theta_theta)**(-1)
-
- ! some useful quantities (see papers for definition of u)
- cosu=Z_theta/dlp
- sinu=-R_theta/dlp
-
- !Jacobian J_r = (dPsi/dr) J_psi = (dPsi/dr) / [(nabla phi x nabla psi)* nabla theta]
- ! = R * (dR/drho dZ/dtheta - dR/dtheta dZ/drho) = R dlp / |nabla r|
- J_r=R*(R_rho*Z_theta-R_theta*Z_rho)
-
- !From definition of q = 1/(2 pi) int (B nabla phi) / (B nabla theta) dtheta:
- !dPsi/dr = sign_Bt sign_Ip / (2 pi q) int F / R^2 J_r dtheta
- ! = F / (2 pi |q|) int J_r/R^2 dtheta
- tmp_arr=J_r/R**2
- drPsi=sign_Ip_CW*F/(2.*pi*n_pol_ext*q0)*sum(tmp_arr)*2*pi*n_pol_ext/np !dlp_int(tmp_arr,1.0)
-
- !Poloidal field (Bp = Bvec * nabla l)
- Bp=sign_Ip_CW * drPsi / J_r * dlp
-
- !toroidal field
- Bphi=F/R
-
- !total modulus of Bfield
- B=sqrt(Bphi**2+Bp**2)
-
- bMaxShift = .false.
- if (thetaShift==0.0.and.trim(magn_geometry).ne.'miller_general') then
- do i = 2,np-1
- if (B(iBmax)<B(i)) then
- iBmax = i
- end if
- enddo
- if (iBmax/=1) then
- bMaxShift = .true.
- thetaShift = theta(iBmax)-theta(1)
- end if
- end if
-
- !definition of radial coordinate! dx_drho=1 --> x = r
- dx_drho=1. !drPsi/Psi0*Lnorm*q0
- if (mype==0.and.print_ini_msg) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho
- dxPsi=drPsi/dx_drho
- geom%C_y=dxPsi*sign_Ip_CW
- geom%Cyq0_x0=geom%C_y(gpdisc%pi1gl)*q0/rho
- geom%C_xy=abs(B0*dxPsi/geom%C_y)
-
- if (mype==0.and.print_ini_msg) then
- write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") &
- "Setting C_xy = ",geom%C_xy(gpdisc%pi1gl),' C_y = ', geom%C_y(gpdisc%pi1gl)," C_x' = ", 1./dxPsi
- write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi
- write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi
- if (thetaShift.ne.0.0) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift
- endif
-
-
- !--------shear is expected to be defined as rho/q*dq/drho--------!
- dq_dx=shat*q0/rho/dx_drho
- dq_dpsi=dq_dx/dxPsi
- pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi
-
- !neg. dpdx normalized to magnetic pressure for pressure term
- geom%dpdx_pm_geom=amhd/q0**2/R0/dx_drho
-
- !first coefficient of psi in varrho expansion
- psi1 = R*Bp*sign_Ip_CW
-
- !integrals for ffprime evaluation
- do i=1,np
- tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
- D0(i)=-F*dlp_int_ind(tmp_arr,dlp,i)
- tmp_arr=B**2*R/psi1**3
- D1(i)=-dlp_int_ind(tmp_arr,dlp,i)/F
- tmp_arr=mu_0*R/psi1**3
- D2(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
- tmp_arr=1./(R*psi1)
- D3(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
- enddo
- tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
- D0_full=-F*dlp_int(tmp_arr,dlp)
- tmp_arr=B**2*R/psi1**3
- D1_full=-dlp_int(tmp_arr,dlp)/F
- tmp_arr=mu_0*R/psi1**3
- D2_full=-dlp_int(tmp_arr,dlp)*F
- tmp_arr=1./(R*psi1)
- D3_full=-dlp_int(tmp_arr,dlp)*F
- D0_mid=D0(np/2+1)
- D1_mid=D1(np/2+1)
- D2_mid=D2(np/2+1)
- D3_mid=D3(np/2+1)
-
- ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*n_pol_ext+D0_full+D2_full*pprime)/D1_full
-
- if (mype==0.and.print_ini_msg) then
- write(*,'(A,ES12.4)') "ffprime = ", ffprime
- endif
- D0=D0-D0_mid
- D1=D1-D1_mid
- D2=D2-D2_mid
- nu=D3-D3_mid
-
- nu1=psi1*(D0+D1*ffprime+D2*pprime)
-
- !straight field line angle defined on equidistant theta grid
- !alpha = phi + nu = - (q chi - phi) => chi = -nu / q
- chi=-nu/q0
-
- !correct small scaling error (<0.5%, due to finite integration resolution)
- chi=chi*(maxval(theta)-minval(theta))/(maxval(chi)-minval(chi))
-
- !new grid equidistant in straight field line angle
- chi_s = linspace(-pi*n_pol_s,pi*n_pol_s-2*pi*n_pol_s/np_s,np_s)
-
- if (sign_Ip_CW.lt.0.0) then !make chi increasing function to not confuse lag3interp
- tmp_reverse = chi(np:1:-1)
- theta_reverse = theta(np:1:-1)
- call lag3interp(theta_reverse,tmp_reverse,np,theta_s,chi_s,np_s)
- theta_s_reverse = theta_s(np_s:1:-1)
- else
- !lag3interp(y_in,x_in,n_in,y_out,x_out,n_out)
- call lag3interp(theta,chi,np,theta_s,chi_s,np_s)
- endif
- dtheta_dchi_s=deriv_fd(theta_s,chi_s,np_s)
-
- !arrays equidistant in straight field line angle
- thAdj_s = theta_s + thetaShift
-
- select case (magn_geometry)
- case ('miller')
- if (zeta/=0.0 .or. s_zeta/=0.0) then
- R_s = R0 + rho*Cos(thAdj_s + d_inv*Sin(thAdj_s))
- Z_s = Z0 + kappa*rho*Sin(thAdj_s + zeta*Sin(2*thAdj_s))
-
- R_theta_s = -(dtheta_dchi_s*rho*(1 + d_inv*Cos(thAdj_s))*Sin(thAdj_s + d_inv*Sin(thAdj_s)))
- Z_theta_s = dtheta_dchi_s*kappa*rho*(1 + 2*zeta*Cos(2*thAdj_s))*Cos(thAdj_s + zeta*Sin(2*thAdj_s))
- else
- Rcirc_s = rho*Cos(thAdj_s - thetad + thetak)
- Zcirc_s = rho*Sin(thAdj_s - thetad + thetak)
- Relong_s = Rcirc_s
- Zelong_s = Zcirc_s + (-1 + kappa)*rho*Sin(thAdj_s - thetad + thetak)
- RelongTilt_s = Relong_s*Cos(thetad - thetak) - Zelong_s*Sin(thetad - thetak)
- ZelongTilt_s = Zelong_s*Cos(thetad - thetak) + Relong_s*Sin(thetad - thetak)
- Rtri_s = RelongTilt_s - rho*Cos(thAdj_s) + rho*Cos(thAdj_s + delta*Sin(thAdj_s))
- Ztri_s = ZelongTilt_s
- RtriTilt_s = Rtri_s*Cos(thetad) + Ztri_s*Sin(thetad)
- ZtriTilt_s = Ztri_s*Cos(thetad) - Rtri_s*Sin(thetad)
- R_s = R0 + RtriTilt_s
- Z_s = Z0 + ZtriTilt_s
-
- dtRcirc_s = -(rho*Sin(thAdj_s - thetad + thetak))
- dtZcirc_s = rho*Cos(thAdj_s - thetad + thetak)
- dtRelong_s = dtRcirc_s
- dtZelong_s = dtZcirc_s + (-1 + kappa)*rho*Cos(thAdj_s - thetad + thetak)
- dtRelongTilt_s = dtRelong_s*Cos(thetad - thetak) - dtZelong_s*Sin(thetad - thetak)
- dtZelongTilt_s = dtZelong_s*Cos(thetad - thetak) + dtRelong_s*Sin(thetad - thetak)
- dtRtri_s = dtRelongTilt_s + rho*Sin(thAdj_s) &
- - rho*(1 + delta*Cos(thAdj_s))*Sin(thAdj_s + delta*Sin(thAdj_s))
- dtZtri_s = dtZelongTilt_s
- dtRtriTilt_s = dtRtri_s*Cos(thetad) + dtZtri_s*Sin(thetad)
- dtZtriTilt_s = dtZtri_s*Cos(thetad) - dtRtri_s*Sin(thetad)
- R_theta_s = dtheta_dchi_s*dtRtriTilt_s
- Z_theta_s = dtheta_dchi_s*dtZtriTilt_s
- endif
- case default
- write (*,*) "ERROR: invalid analytic geometry specification"
- end select
-
- if (sign_Ip_CW.lt.0.0) then
- call lag3interp(nu1,theta,np,tmp_arr_s,theta_s_reverse,np_s)
- nu1_s = tmp_arr_s(np_s:1:-1)
- call lag3interp(Bp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
- Bp_s = tmp_arr_s(np_s:1:-1)
- call lag3interp(dlp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
- dlp_s = tmp_arr_s(np_s:1:-1)
- call lag3interp(Rc,theta,np,tmp_arr_s,theta_s_reverse,np_s)
- Rc_s = tmp_arr_s(np_s:1:-1)
- else
- call lag3interp(nu1,theta,np,nu1_s,theta_s,np_s)
- call lag3interp(Bp,theta,np,Bp_s,theta_s,np_s)
- call lag3interp(dlp,theta,np,dlp_s,theta_s,np_s)
- call lag3interp(Rc,theta,np,Rc_s,theta_s,np_s)
- endif
-
- psi1_s = R_s*Bp_s*sign_Ip_CW
-
- dBp_dchi_s=deriv_fd(Bp_s,chi_s,np_s)
-
- Bphi_s=F/R_s
- B_s=sqrt(Bphi_s**2+Bp_s**2)
- cosu_s=Z_theta_s/dlp_s/dtheta_dchi_s
- sinu_s=-R_theta_s/dlp_s/dtheta_dchi_s
-
- !radial derivative of straight field line angle
- dchidx_s=-(nu1_s/psi1_s*dxPsi+chi_s*dq_dx)/q0
-
- !Bfield derivatives in Mercier-Luc coordinates (varrho,l,phi)
- dB_drho_s=-1./B_s*(F**2/R_s**3*cosu_s+Bp_s**2/Rc_s+mu_0*psi1_s*pprime)
- dB_dl_s=1./B_s*(Bp_s*dBp_dchi_s/dtheta_dchi_s/dlp_s+F**2/R_s**3*sinu_s)
-
- dnu_drho_s=nu1_s
- dnu_dl_s=-F/(R_s*psi1_s)
- grad_nu_s=sqrt(dnu_drho_s**2+dnu_dl_s**2)
-
- !contravariant metric coefficients (varrho,l,phi)->(x,y,z)
- gxx=(psi1_s/dxPsi)**2
- gxy=-psi1_s/dxPsi*geom%C_y(gpdisc%pi1gl)*sign_Ip_CW*nu1_s
- gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0
- gyy=geom%C_y(gpdisc%pi1gl)**2*(grad_nu_s**2+1/R_s**2)
- gyz=sign_Ip_CW*geom%C_y(gpdisc%pi1gl)/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s)
- gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2)
-
- jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2)
-
- !covariant metric coefficients
- g_xx=jacobian**2*(gyy*gzz-gyz**2)
- g_xy=jacobian**2*(gxz*gyz-gxy*gzz)
- g_xz=jacobian**2*(gxy*gyz-gxz*gyy)
- g_yy=jacobian**2*(gxx*gzz-gxz**2)
- g_yz=jacobian**2*(gxz*gxy-gxx*gyz)
- g_zz=jacobian**2*(gxx*gyy-gxy**2)
-
- !Bfield derivatives
- !dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B
- dBdx=jacobian*geom%C_y(gpdisc%pi1gl)/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s)
- dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s)
-
- !curvature terms (these are just local and will be recalculated in geometry.F90)
- K_x = (0.-g_yz/g_zz*dBdz)
- K_y = (dBdx-g_xz/g_zz*dBdz)
-
- !(R,Z) derivatives for visualization
- dxdR_s = dx_drho/drPsi*psi1_s*cosu_s
- dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s
-
- if (edge_opt==0.0) then
- !gene z-grid
- chi_out=linspace(-pi*n_pol,pi*n_pol-2*pi*n_pol/nz0,nz0)
- else
- !new parallel coordinate chi_out==zprime
- !see also tracer_aux.F90
- if (n_pol>1) STOP "ERROR: n_pol>1 has not been implemented for edge_opt=\=0.0"
- do k=0,nz0-1
- chi_out(k)=sinh((-pi+k*2.*pi/nz0)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
- enddo
- !transform metrics according to chain rule
- do k=1,np_s
- dzprimedz=get_dzprimedz(chi_s(k))
- gxz(k)=gxz(k)*dzprimedz
- gyz(k)=gyz(k)*dzprimedz
- gzz(k)=gzz(k)*dzprimedz**2
- jacobian(k)=jacobian(k)/dzprimedz
- dBdz(k)=dBdz(k)/dzprimedz
- enddo
- endif !edge_opt
-
- !interpolate down to GENE z-grid
- call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,nz0)
- call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,nz0)
- call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,nz0)
- call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,nz0)
- call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,nz0)
- call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,nz0)
- call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,nz0)
- call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,nz0)
- call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,nz0)
- call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,nz0)
- call lag3interp(R_s,chi_s,np_s,R_out,chi_out,nz0)
- call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,nz0)
- call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,nz0)
- call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,nz0)
-
- !select local k range
- do i=gpdisc%pi1gl,gpdisc%pi2gl
- if (gdisc%yx_order) then
- geom%gii(i,lk1:lk2)=gyy_out(lk1:lk2)
- geom%gjj(i,lk1:lk2)=gxx_out(lk1:lk2)
- geom%giz(i,lk1:lk2)=gyz_out(lk1:lk2)
- geom%gjz(i,lk1:lk2)=gxz_out(lk1:lk2)
- geom%dBdi(i,lk1:lk2)=0.
- geom%dBdj(i,lk1:lk2)=dBdx_out(lk1:lk2)
- else
- geom%gii(i,lk1:lk2)=gxx_out(lk1:lk2)
- geom%gjj(i,lk1:lk2)=gyy_out(lk1:lk2)
- geom%giz(i,lk1:lk2)=gxz_out(lk1:lk2)
- geom%gjz(i,lk1:lk2)=gyz_out(lk1:lk2)
- geom%dBdi(i,lk1:lk2)=dBdx_out(lk1:lk2)
- geom%dBdj(i,lk1:lk2)=0.
- endif
- geom%gij(i,lk1:lk2)=gxy_out(lk1:lk2)
- geom%gzz(i,lk1:lk2)=gzz_out(lk1:lk2)
- geom%Bfield(i,lk1:lk2)=Bfield_out(lk1:lk2)
- geom%jacobian(i,lk1:lk2)=jacobian_out(lk1:lk2)
- geom%dBdz(i,lk1:lk2)=dBdz_out(lk1:lk2)
- if (Lref.gt.0.) then
- geom%R(i,lk1:lk2)=R_out(lk1:lk2)*Lref
- geom%Z(i,lk1:lk2)=Z_out(lk1:lk2)*Lref
- else
- geom%R(i,lk1:lk2)=R_out(lk1:lk2)
- geom%Z(i,lk1:lk2)=Z_out(lk1:lk2)
- endif
- geom%R_hat(i,lk1:lk2)=R_out(lk1:lk2)
- geom%Z_hat(i,lk1:lk2)=Z_out(lk1:lk2)
- geom%dxdR(i,lk1:lk2)= dxdR_out(lk1:lk2)
- geom%dxdZ(i,lk1:lk2)= dxdZ_out(lk1:lk2)
- enddo
-
- geom%q_prof=q0
- geom%dqdx_prof=shat*q0/gcoor%x0
-
- contains
-
- !> Generate an equidistant array from min to max with n points
- function linspace(min,max,n) result(out)
- real:: min, max
- integer:: n
- real, dimension(n):: out
-
- do i=1,n
- out(i)=min+(i-1)*(max-min)/(n-1)
- enddo
- end function linspace
-
- !> Weighted average
- real function average(var,weight)
- real, dimension(np):: var, weight
- average=sum(var*weight)/sum(weight)
- end function average
-
- !> full theta integral with weight function dlp
- real function dlp_int(var,dlp)
- real, dimension(np):: var, dlp
- dlp_int=sum(var*dlp)*2*pi*n_pol_ext/np
- end function dlp_int
-
- !> theta integral with weight function dlp, up to index 'ind'
- real function dlp_int_ind(var,dlp,ind)
- real, dimension(np):: var, dlp
- integer:: ind
-
- dlp_int_ind=0.
- if (ind.gt.1) then
- dlp_int_ind=dlp_int_ind+var(1)*dlp(1)*pi*n_pol_ext/np
- dlp_int_ind=dlp_int_ind+(sum(var(2:ind-1)*dlp(2:ind-1)))*2*pi*n_pol_ext/np
- dlp_int_ind=dlp_int_ind+var(ind)*dlp(ind)*pi*n_pol_ext/np
- endif
- end function dlp_int_ind
-
- !> 1st derivative with 2nd order finite differences
- function deriv_fd(y,x,n) result(out)
- integer, intent(in) :: n
- real, dimension(n):: x,y,out,dx
-
- !call lag3deriv(y,x,n,out,x,n)
-
- out=0.
- do i=2,n-1
- out(i)=out(i)-y(i-1)/2
- out(i)=out(i)+y(i+1)/2
- enddo
- out(1)=y(2)-y(1)
- out(n)=y(n)-y(n-1)
- dx=x(2)-x(1)
- out=out/dx
-
- end function deriv_fd
-
-
- end subroutine get_miller
-
-
-END MODULE miller_mod
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
new file mode 100644
index 00000000..f75a075e
--- /dev/null
+++ b/src/miller_mod.F90
@@ -0,0 +1,645 @@
+!! This source has been adapted from GENE https://genecode.org/ !!
+!>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998)
+!>and [J. Candy, PPCF 51, 105009 (2009)]
+MODULE miller
+ USE prec_const
+ USE basic
+ ! use coordinates,only: gcoor, get_dzprimedz
+ USE grid
+ ! use discretization
+ USE lagrange_interpolation
+ ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol
+ USE model
+
+ implicit none
+ public:: get_miller, set_miller_parameters
+ public:: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
+ public:: thetaShift
+ public:: mMode, nMode
+ public:: thetak, thetad
+ public:: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
+ public:: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
+ public:: cN_m, sN_m, cNdr_m, sNdr_m
+
+ private
+
+ real(dp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
+
+ INTEGER :: mMode, nMode
+ real(dp) :: thetaShift
+ real(dp) :: thetak, thetad
+ real(dp) :: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
+ real(dp) :: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
+
+ INTEGER, PARAMETER :: IND_M=32
+ real(dp), DIMENSION(0:IND_M-1) :: cN_m, sN_m, cNdr_m, sNdr_m
+
+CONTAINS
+
+ !>Set defaults for miller parameters
+ subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_)
+ real(dp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_
+ rho = -1.0
+ kappa = kappa_
+ s_kappa = s_kappa_
+ delta = delta_
+ s_delta = s_delta_
+ zeta = zeta_
+ s_zeta = s_zeta_
+ drR = 0.0
+ drZ = 0.0
+
+ thetak = 0.0
+ thetad = 0.0
+
+ aSurf = 0.54
+ Delta2 = 1.0
+ Delta3 = 1.0
+ theta2 = 0.0
+ theta3 = 0.0
+ Raxis = 1.0
+ Zaxis = 0.0
+
+ mMode = 2
+ nMode = 3
+ Deltam = 1.0
+ Deltan = 1.0
+ s_Deltam = 0.0
+ s_Deltan = 0.0
+ thetam = 0.0
+ thetan = 0.0
+
+ cN_m = 0.0
+ sN_m = 0.0
+ cNdr_m = 0.0
+ sNdr_m = 0.0
+ end subroutine set_miller_parameters
+
+ !>Get Miller metric, magnetic field, jacobian etc.
+ subroutine get_miller(trpeps,major_R,major_Z,q0,shat,amhd,edge_opt,&
+ C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,&
+ Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_)
+ !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
+ real(dp), INTENT(INOUT) :: major_R ! major radius
+ real(dp), INTENT(INOUT) :: major_Z ! major Z
+ real(dp), INTENT(INOUT) :: q0 ! safetyfactor
+ real(dp), INTENT(INOUT) :: shat ! safetyfactor
+ real(dp), INTENT(INOUT) :: amhd ! alpha mhd
+ real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd
+ real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad.
+ real(dp), INTENT(INOUT) :: C_y, C_xy
+
+ real(dp), dimension(izgs:izge,0:1), INTENT(INOUT) :: &
+ gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
+ dBdx_,dBdy_,Bfield_,jacobian_,&
+ dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_
+ ! No parameter in gyacomo yet
+ integer :: ikxgs =1 ! the left ghost gpdisc%pi1gl
+ real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
+ real(dp) :: sign_Bt_CW=1 ! current sign (only normal current)
+ !!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ integer:: np, np_s, Npol_ext, Npol_s
+
+ real(dp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi
+ real(dp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt
+ real(dp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt
+ real(dp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
+ real(dp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
+ real(dp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
+ real(dp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
+ real(dp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
+ real(dp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
+ real(dp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
+
+ real(dp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse
+ real(dp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s
+ real(dp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s
+ real(dp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s
+ real(dp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s, dtrShape_s
+ real(dp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s
+ real(dp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s
+ real(dp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
+ real(dp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
+
+ real(dp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
+ real(dp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out
+ real(dp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
+ !real(dp) :: Lnorm, Psi0 ! currently module-wide defined anyway
+ real(dp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
+ real(dp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
+ logical:: bMaxShift
+ integer:: i, k, iBmax
+
+ Npol_ext = Npol+2
+ Npol_s = Npol+1
+ np = 500*Npol_ext
+ np_s = 500*Npol_s
+
+ rho = trpeps*major_R
+ if (rho.le.0.0) stop 'flux surface radius not defined'
+ trpeps = rho/major_R
+
+ q0 = sign_Ip_CW * sign_Bt_CW * abs(q0)
+
+ R0=major_R
+ B0=1.0*sign_Bt_CW
+ F=R0*B0
+ Z0=major_Z
+ pi = acos(-1.0)
+ mu_0=4.0*pi
+
+ theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._dp*pi*Npol_ext/np,np)
+ d_inv=asin(delta)
+
+ thetaShift = 0.0
+ iBmax = 1
+
+ !flux surface parametrization
+ thAdj = theta + thetaShift
+ if (zeta/=0.0 .or. s_zeta/=0.0) then
+ R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
+ Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
+ Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
+ + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
+ else
+ Rcirc = rho*Cos(thAdj - thetad + thetak)
+ Zcirc = rho*Sin(thAdj - thetad + thetak)
+ Relong = Rcirc
+ Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
+ ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
+ Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
+ Ztri = ZelongTilt
+ RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
+ ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
+ R = R0 + RtriTilt
+ Z = Z0 + ZtriTilt
+
+ drRcirc = Cos(thAdj - thetad + thetak)
+ drZcirc = Sin(thAdj - thetad + thetak)
+ drRelong = drRcirc
+ drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
+ drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
+ drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
+ drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
+ - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ drZtri = drZelongTilt
+ drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
+ drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
+ R_rho = drR + drRtriTilt
+ Z_rho = drZ + drZtriTilt
+
+ dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtZcirc = rho*Cos(thAdj - thetad + thetak)
+ dtRelong = dtRcirc
+ dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
+ dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
+ dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
+ dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
+ dtZtri = dtZelongTilt
+ dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
+ dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
+ R_theta = dtRtriTilt
+ Z_theta = dtZtriTilt
+
+ dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
+ dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtdtRelong = dtdtRcirc
+ dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
+ dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
+ dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
+ + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ dtdtZtri = dtdtZelongTilt
+ dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
+ dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
+ R_theta_theta = dtdtRtriTilt
+ Z_theta_theta = dtdtZtriTilt
+ endif
+
+ !dl/dtheta
+ dlp=(R_theta**2+Z_theta**2)**0.5
+
+ !curvature radius
+ Rc=dlp**3*(R_theta*Z_theta_theta-Z_theta*R_theta_theta)**(-1)
+
+ ! some useful quantities (see papers for definition of u)
+ cosu=Z_theta/dlp
+ sinu=-R_theta/dlp
+
+ !Jacobian J_r = (dPsi/dr) J_psi = (dPsi/dr) / [(nabla fz x nabla psi)* nabla theta]
+ ! = R * (dR/drho dZ/dtheta - dR/dtheta dZ/drho) = R dlp / |nabla r|
+ J_r=R*(R_rho*Z_theta-R_theta*Z_rho)
+
+ !From definition of q = 1/(2 pi) int (B nabla fz) / (B nabla theta) dtheta:
+ !dPsi/dr = sign_Bt sign_Ip / (2 pi q) int F / R^2 J_r dtheta
+ ! = F / (2 pi |q|) int J_r/R^2 dtheta
+ tmp_arr=J_r/R**2
+ drPsi=sign_Ip_CW*F/(2.*pi*Npol_ext*q0)*sum(tmp_arr)*2*pi*Npol_ext/np !dlp_int(tmp_arr,1.0)
+
+ !Poloidal field (Bp = Bvec * nabla l)
+ Bp=sign_Ip_CW * drPsi / J_r * dlp
+
+ !toroidal field
+ Bphi=F/R
+
+ !total modulus of Bfield
+ B=sqrt(Bphi**2+Bp**2)
+
+ bMaxShift = .false.
+ ! if (thetaShift==0.0.and.trim(magn_geometry).ne.'miller_general') then
+ if (thetaShift==0.0) then
+ do i = 2,np-1
+ if (B(iBmax)<B(i)) then
+ iBmax = i
+ end if
+ enddo
+ if (iBmax/=1) then
+ bMaxShift = .true.
+ thetaShift = theta(iBmax)-theta(1)
+ end if
+ end if
+
+ !definition of radial coordinate! dx_drho=1 --> x = r
+ dx_drho=1. !drPsi/Psi0*Lnorm*q0
+ if (my_id==0) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho
+ dxPsi = drPsi/dx_drho
+ C_y = dxPsi*sign_Ip_CW
+ C_xy = abs(B0*dxPsi/C_y)
+
+ if (my_id==0) then
+ write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") &
+ "Setting C_xy = ",C_xy,' C_y = ', C_y," C_x' = ", 1./dxPsi
+ write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi
+ write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi
+ if (thetaShift.ne.0.0) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift
+ endif
+
+
+ !--------shear is expected to be defined as rho/q*dq/drho--------!
+ dq_dx=shat*q0/rho/dx_drho
+ dq_dpsi=dq_dx/dxPsi
+ pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi
+
+ !neg. dpdx normalized to magnetic pressure for pressure term
+ dpdx_pm_geom=amhd/q0**2/R0/dx_drho
+
+ !first coefficient of psi in varrho expansion
+ psi1 = R*Bp*sign_Ip_CW
+
+ !integrals for ffprime evaluation
+ do i=1,np
+ tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
+ D0(i)=-F*dlp_int_ind(tmp_arr,dlp,i)
+ tmp_arr=B**2*R/psi1**3
+ D1(i)=-dlp_int_ind(tmp_arr,dlp,i)/F
+ tmp_arr=mu_0*R/psi1**3
+ D2(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
+ tmp_arr=1./(R*psi1)
+ D3(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
+ enddo
+ tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
+ D0_full=-F*dlp_int(tmp_arr,dlp)
+ tmp_arr=B**2*R/psi1**3
+ D1_full=-dlp_int(tmp_arr,dlp)/F
+ tmp_arr=mu_0*R/psi1**3
+ D2_full=-dlp_int(tmp_arr,dlp)*F
+ tmp_arr=1./(R*psi1)
+ D3_full=-dlp_int(tmp_arr,dlp)*F
+ D0_mid=D0(np/2+1)
+ D1_mid=D1(np/2+1)
+ D2_mid=D2(np/2+1)
+ D3_mid=D3(np/2+1)
+
+ ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full
+
+ if (my_id==0) then
+ write(*,'(A,ES12.4)') "ffprime = ", ffprime
+ endif
+ D0=D0-D0_mid
+ D1=D1-D1_mid
+ D2=D2-D2_mid
+ nu=D3-D3_mid
+
+ nu1=psi1*(D0+D1*ffprime+D2*pprime)
+
+ !straight field line angle defined on equidistant theta grid
+ !alpha = fz + nu = - (q chi - fz) => chi = -nu / q
+ chi=-nu/q0
+
+ !correct small scaling error (<0.5%, due to finite integration resolution)
+ chi=chi*(maxval(theta)-minval(theta))/(maxval(chi)-minval(chi))
+
+ !new grid equidistant in straight field line angle
+ chi_s = linspace(-pi*Npol_s,pi*Npol_s-2*pi*Npol_s/np_s,np_s)
+
+ if (sign_Ip_CW.lt.0.0) then !make chi increasing function to not confuse lag3interp
+ tmp_reverse = chi(np:1:-1)
+ theta_reverse = theta(np:1:-1)
+ call lag3interp(theta_reverse,tmp_reverse,np,theta_s,chi_s,np_s)
+ theta_s_reverse = theta_s(np_s:1:-1)
+ else
+ !lag3interp(y_in,x_in,n_in,y_out,x_out,n_out)
+ call lag3interp(theta,chi,np,theta_s,chi_s,np_s)
+ endif
+ dtheta_dchi_s=deriv_fd(theta_s,chi_s,np_s)
+
+ !arrays equidistant in straight field line angle
+ thAdj_s = theta_s + thetaShift
+
+ if (zeta/=0.0 .or. s_zeta/=0.0) then
+ R_s = R0 + rho*Cos(thAdj_s + d_inv*Sin(thAdj_s))
+ Z_s = Z0 + kappa*rho*Sin(thAdj_s + zeta*Sin(2*thAdj_s))
+
+ R_theta_s = -(dtheta_dchi_s*rho*(1 + d_inv*Cos(thAdj_s))*Sin(thAdj_s + d_inv*Sin(thAdj_s)))
+ Z_theta_s = dtheta_dchi_s*kappa*rho*(1 + 2*zeta*Cos(2*thAdj_s))*Cos(thAdj_s + zeta*Sin(2*thAdj_s))
+ else
+ Rcirc_s = rho*Cos(thAdj_s - thetad + thetak)
+ Zcirc_s = rho*Sin(thAdj_s - thetad + thetak)
+ Relong_s = Rcirc_s
+ Zelong_s = Zcirc_s + (-1 + kappa)*rho*Sin(thAdj_s - thetad + thetak)
+ RelongTilt_s = Relong_s*Cos(thetad - thetak) - Zelong_s*Sin(thetad - thetak)
+ ZelongTilt_s = Zelong_s*Cos(thetad - thetak) + Relong_s*Sin(thetad - thetak)
+ Rtri_s = RelongTilt_s - rho*Cos(thAdj_s) + rho*Cos(thAdj_s + delta*Sin(thAdj_s))
+ Ztri_s = ZelongTilt_s
+ RtriTilt_s = Rtri_s*Cos(thetad) + Ztri_s*Sin(thetad)
+ ZtriTilt_s = Ztri_s*Cos(thetad) - Rtri_s*Sin(thetad)
+ R_s = R0 + RtriTilt_s
+ Z_s = Z0 + ZtriTilt_s
+
+ dtRcirc_s = -(rho*Sin(thAdj_s - thetad + thetak))
+ dtZcirc_s = rho*Cos(thAdj_s - thetad + thetak)
+ dtRelong_s = dtRcirc_s
+ dtZelong_s = dtZcirc_s + (-1 + kappa)*rho*Cos(thAdj_s - thetad + thetak)
+ dtRelongTilt_s = dtRelong_s*Cos(thetad - thetak) - dtZelong_s*Sin(thetad - thetak)
+ dtZelongTilt_s = dtZelong_s*Cos(thetad - thetak) + dtRelong_s*Sin(thetad - thetak)
+ dtRtri_s = dtRelongTilt_s + rho*Sin(thAdj_s) &
+ - rho*(1 + delta*Cos(thAdj_s))*Sin(thAdj_s + delta*Sin(thAdj_s))
+ dtZtri_s = dtZelongTilt_s
+ dtRtriTilt_s = dtRtri_s*Cos(thetad) + dtZtri_s*Sin(thetad)
+ dtZtriTilt_s = dtZtri_s*Cos(thetad) - dtRtri_s*Sin(thetad)
+ R_theta_s = dtheta_dchi_s*dtRtriTilt_s
+ Z_theta_s = dtheta_dchi_s*dtZtriTilt_s
+ endif
+ if (sign_Ip_CW.lt.0.0) then
+ call lag3interp(nu1,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ nu1_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(Bp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ Bp_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(dlp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ dlp_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(Rc,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ Rc_s = tmp_arr_s(np_s:1:-1)
+ else
+ call lag3interp(nu1,theta,np,nu1_s,theta_s,np_s)
+ call lag3interp(Bp,theta,np,Bp_s,theta_s,np_s)
+ call lag3interp(dlp,theta,np,dlp_s,theta_s,np_s)
+ call lag3interp(Rc,theta,np,Rc_s,theta_s,np_s)
+ endif
+
+ psi1_s = R_s*Bp_s*sign_Ip_CW
+
+ dBp_dchi_s=deriv_fd(Bp_s,chi_s,np_s)
+
+ Bphi_s=F/R_s
+ B_s=sqrt(Bphi_s**2+Bp_s**2)
+ cosu_s=Z_theta_s/dlp_s/dtheta_dchi_s
+ sinu_s=-R_theta_s/dlp_s/dtheta_dchi_s
+
+ !radial derivative of straight field line angle
+ dchidx_s=-(nu1_s/psi1_s*dxPsi+chi_s*dq_dx)/q0
+
+ !Bfield derivatives in Mercier-Luc coordinates (varrho,l,fz)
+ dB_drho_s=-1./B_s*(F**2/R_s**3*cosu_s+Bp_s**2/Rc_s+mu_0*psi1_s*pprime)
+ dB_dl_s=1./B_s*(Bp_s*dBp_dchi_s/dtheta_dchi_s/dlp_s+F**2/R_s**3*sinu_s)
+
+ dnu_drho_s=nu1_s
+ dnu_dl_s=-F/(R_s*psi1_s)
+ grad_nu_s=sqrt(dnu_drho_s**2+dnu_dl_s**2)
+
+ !contravariant metric coefficients (varrho,l,fz)->(x,y,z)
+ gxx=(psi1_s/dxPsi)**2
+ gxy=-psi1_s/dxPsi*C_y*sign_Ip_CW*nu1_s
+ gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0
+ gyy=C_y**2*(grad_nu_s**2+1/R_s**2)
+ gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s)
+ gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2)
+
+ jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2)
+
+ !covariant metric coefficients
+ g_xx=jacobian**2*(gyy*gzz-gyz**2)
+ g_xy=jacobian**2*(gxz*gyz-gxy*gzz)
+ g_xz=jacobian**2*(gxy*gyz-gxz*gyy)
+ g_yy=jacobian**2*(gxx*gzz-gxz**2)
+ g_yz=jacobian**2*(gxz*gxy-gxx*gyz)
+ g_zz=jacobian**2*(gxx*gyy-gxy**2)
+
+ !Bfield derivatives
+ !dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B
+ dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s)
+ dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s)
+
+ !curvature terms (these are just local and will be recalculated in geometry.F90)
+ K_x = (0.-g_yz/g_zz*dBdz)
+ K_y = (dBdx-g_xz/g_zz*dBdz)
+
+ !(R,Z) derivatives for visualization
+ dxdR_s = dx_drho/drPsi*psi1_s*cosu_s
+ dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s
+
+ if (edge_opt==0.0) then
+ !gene z-grid
+ chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/Nz,Nz)
+ else
+ !new parallel coordinate chi_out==zprime
+ !see also tracer_aux.F90
+ if (Npol>1) STOP "ERROR: Npol>1 has not been implemented for edge_opt=\=0.0"
+ do k=izs,ize
+ chi_out(k)=sinh((-pi+k*2.*pi/Nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
+ enddo
+ !transform metrics according to chain rule
+ do k=1,np_s
+ !>dz'/dz conversion for edge_opt as function of z
+ if (edge_opt.gt.0) then
+ dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1))/&
+ sqrt((edge_opt*chi_s(k))**2+1)
+ else
+ dzprimedz = 1.0
+ endif
+ gxz(k)=gxz(k)*dzprimedz
+ gyz(k)=gyz(k)*dzprimedz
+ gzz(k)=gzz(k)*dzprimedz**2
+ jacobian(k)=jacobian(k)/dzprimedz
+ dBdz(k)=dBdz(k)/dzprimedz
+ enddo
+ endif !edge_opt
+
+ !interpolate down to GENE z-grid
+ call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,Nz)
+ call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,Nz)
+ call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,Nz)
+ call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,Nz)
+ call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,Nz)
+ call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,Nz)
+ call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,Nz)
+ call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,Nz)
+ call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,Nz)
+ call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,Nz)
+ call lag3interp(R_s,chi_s,np_s,R_out,chi_out,Nz)
+ call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,Nz)
+ call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,Nz)
+ call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz)
+ ! Fill the geom arrays with the results
+ do eo=0,1
+ gxx_(izs:ize,eo) =gxx_out(izs:ize)
+ gyy_(izs:ize,eo) =gyy_out(izs:ize)
+ gxz_(izs:ize,eo) =gxz_out(izs:ize)
+ gyz_(izs:ize,eo) =gyz_out(izs:ize)
+ dBdx_(izs:ize,eo) =dBdx_out(izs:ize)
+ dBdy_(izs:ize,eo) =0.
+ gxy_(izs:ize,eo) =gxy_out(izs:ize)
+ gzz_(izs:ize,eo) =gzz_out(izs:ize)
+ Bfield_(izs:ize,eo) =Bfield_out(izs:ize)
+ jacobian_(izs:ize,eo) =jacobian_out(izs:ize)
+ dBdz_(izs:ize,eo) =dBdz_out(izs:ize)
+ R_hat_(izs:ize,eo) =R_out(izs:ize)
+ Z_hat_(izs:ize,eo) =Z_out(izs:ize)
+ dxdR_(izs:ize,eo) = dxdR_out(izs:ize)
+ dxdZ_(izs:ize,eo) = dxdZ_out(izs:ize)
+
+ !! UPDATE GHOSTS VALUES (since the miller function in GENE does not)
+ CALL update_ghosts_z(gxx_(:,eo))
+ CALL update_ghosts_z(gyy_(:,eo))
+ CALL update_ghosts_z(gxz_(:,eo))
+ CALL update_ghosts_z(dBdx_(:,eo))
+ CALL update_ghosts_z(dBdy_(:,eo))
+ CALL update_ghosts_z(gxy_(:,eo))
+ CALL update_ghosts_z(gzz_(:,eo))
+ CALL update_ghosts_z(Bfield_(:,eo))
+ CALL update_ghosts_z(jacobian_(:,eo))
+ CALL update_ghosts_z(dBdz_(:,eo))
+ CALL update_ghosts_z(R_hat_(:,eo))
+ CALL update_ghosts_z(Z_hat_(:,eo))
+ CALL update_ghosts_z(dxdR_(:,eo))
+ CALL update_ghosts_z(dxdZ_(:,eo))
+ enddo
+
+ contains
+
+
+ SUBROUTINE update_ghosts_z(fz_)
+ IMPLICIT NONE
+ ! INTEGER, INTENT(IN) :: nztot_
+ REAL(dp), DIMENSION(izgs:izge), INTENT(INOUT) :: fz_
+ REAL(dp), DIMENSION(-2:2) :: buff
+ INTEGER :: status(MPI_STATUS_SIZE), count
+
+ IF(Nz .GT. 1) THEN
+ IF (num_procs_z .GT. 1) THEN
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ count = 1 ! one point to exchange
+ !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(fz_(ize), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last
+ buff(-1), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv(fz_(ize-1), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last
+ buff(-2), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-2
+ comm0, status, ierr)
+
+ !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(fz_(izs), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first
+ buff(+1), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv(fz_(izs+1), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first
+ buff(+2), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+2
+ comm0, status, ierr)
+ ELSE
+ buff(-1) = fz_(ize )
+ buff(-2) = fz_(ize-1)
+ buff(+1) = fz_(izs )
+ buff(+2) = fz_(izs+1)
+ ENDIF
+ fz_(ize+1) = buff(+1)
+ fz_(ize+2) = buff(+2)
+ fz_(izs-1) = buff(-1)
+ fz_(izs-2) = buff(-2)
+ ENDIF
+ END SUBROUTINE update_ghosts_z
+
+
+ !> Generate an equidistant array from min to max with n points
+ function linspace(min,max,n) result(out)
+ real(dp):: min, max
+ integer:: n
+ real(dp), dimension(n):: out
+
+ do i=1,n
+ out(i)=min+(i-1)*(max-min)/(n-1)
+ enddo
+ end function linspace
+
+ !> Weighted average
+ real(dp) function average(var,weight)
+ real(dp), dimension(np):: var, weight
+ average=sum(var*weight)/sum(weight)
+ end function average
+
+ !> full theta integral with weight function dlp
+ real(dp) function dlp_int(var,dlp)
+ real(dp), dimension(np):: var, dlp
+ dlp_int=sum(var*dlp)*2*pi*Npol_ext/np
+ end function dlp_int
+
+ !> theta integral with weight function dlp, up to index 'ind'
+ real(dp) function dlp_int_ind(var,dlp,ind)
+ real(dp), dimension(np):: var, dlp
+ integer:: ind
+
+ dlp_int_ind=0.
+ if (ind.gt.1) then
+ dlp_int_ind=dlp_int_ind+var(1)*dlp(1)*pi*Npol_ext/np
+ dlp_int_ind=dlp_int_ind+(sum(var(2:ind-1)*dlp(2:ind-1)))*2*pi*Npol_ext/np
+ dlp_int_ind=dlp_int_ind+var(ind)*dlp(ind)*pi*Npol_ext/np
+ endif
+ end function dlp_int_ind
+
+ !> 1st derivative with 2nd order finite differences
+ function deriv_fd(y,x,n) result(out)
+ integer, intent(in) :: n
+ real(dp), dimension(n):: x,y,out,dx
+
+ !call lag3deriv(y,x,n,out,x,n)
+
+ out=0.
+ do i=2,n-1
+ out(i)=out(i)-y(i-1)/2
+ out(i)=out(i)+y(i+1)/2
+ enddo
+ out(1)=y(2)-y(1)
+ out(n)=y(n)-y(n-1)
+ dx=x(2)-x(1)
+ out=out/dx
+
+ end function deriv_fd
+
+
+ end subroutine get_miller
+
+
+END MODULE miller
--
GitLab