module geometry
! computes geometrical quantities
! Adapted from B.J.Frei MOLIX code (2021)
use prec_const
use model
use grid
use array
use fields
use basic
use calculus, ONLY: simpson_rule_z
implicit none
PRIVATE
! Geometry 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
! Curvature
REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
! Metric
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
! Relative magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL
! Relative strength of major radius
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
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
CONTAINS
SUBROUTINE geometry_readinputs
! Read the input parameters
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps
READ(lu_in,geometry)
IF(shear .NE. 0._dp) SHEARED = .true.
END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
INTEGER :: fid
! Allocate arrays
CALL geometry_allocate_mem
!
IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
call eval_1D_geometry
ELSE
SELECT CASE(geom)
CASE('circular')
IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
call eval_circular_geometry
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
CASE('Z-pinch')
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
SHEARED = .FALSE.
CASE DEFAULT
ERROR STOP 'Error stop: geometry not recognized!!'
END SELECT
ENDIF
!
! Evaluate perpendicular wavenumber
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
DO eo = 0,1
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx = ikxs, ikxe
kx = kxarray(ikx)
DO iz = izgs,izge
kparray(iky, ikx, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
ENDDO
ENDDO
ENDDO
ENDDO
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
integrant = Jacobian(izs:ize,0) ! Convert into complex array
CALL simpson_rule_z(integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z, kx, ky, alpha_MHD, Gx, Gy
alpha_MHD = 0._dp
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = shear*z - alpha_MHD*SIN(z)
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
gyz(iz,eo) = 1._dp/eps
gzz(iz,eo) = 0._dp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
hatR(iz,eo) = 1._dp + eps*COS(z)
hatZ(iz,eo) = 1._dp + eps*SIN(z)
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
hatB (iz,eo) = 1._dp / hatR(iz,eo)
hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = eps*SIN(z)/hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Kx
Gy = -COS(z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)* hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_salphaB_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_circular_geometry
! evaluate circular geometry model
! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
implicit none
REAL(dp) :: chi, kx, ky, Gx, Gy
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
chi = zarray(iz,eo) ! = chi
! metric in x,y,z
gxx(iz,eo) = 1._dp
gyy(iz,eo) = 1._dp + (shear*chi)**2 - 2._dp*eps*COS(chi) - 2._dp*shear*chi*eps*SIN(chi)
gxy(iz,eo) = shear*chi - eps*SIN(chi)
gxz(iz,eo) = -SIN(chi)
gyz(iz,eo) = (1._dp - 2._dp*eps*COS(chi) - shear*chi*eps*SIN(chi))/eps
gzz(iz,eo) = (1._dp - 2._dp*eps*COS(chi))/eps**2
dxdR(iz,eo)= COS(chi)
dxdZ(iz,eo)= SIN(chi)
! Relative strengh of radius
hatR(iz,eo) = 1._dp + eps*COS(chi)
hatZ(iz,eo) = 1._dp + eps*SIN(chi)
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = chi
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
hatB (iz,eo) = 1._dp / hatR(iz,eo)
hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(chi) ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = eps * SIN(chi) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Kx
Gy = -COS(chi) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Ky
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_circular_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_circular_geometry_GENE
! evaluate circular geometry model
! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
implicit none
REAL(dp) :: qbar, dxdr_, dpsidr, dqdr, dqbardr, Cy, dCydr_Cy, Cxy, fac2, &
cost, sint, dchidr, dchidt, sign_Ip, B, dBdt, dBdr, &
g11, g22, g12, g33, dBdchi, dBdr_c !GENE variables
REAL(dp) :: X, z, kx, ky, Gamma1, Gamma2, Gamma3, feps
sign_Ip = 1._dp
feps = SQRT(1._dp-eps**2)
qbar = q0 * feps
dxdr_ = 1._dp
dpsidr = eps/qbar
dqdr = q0/eps*shear
dqbardr = feps*q0/eps*(shear - eps**2/feps**2)
Cy = eps/q0 * sign_Ip
dCydr_Cy = 0._dp
Cxy = 1._dp/feps
fac2 = dCydr_Cy + dqdr/q0
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
z = zarray(iz,eo)
cost = (COS(z)-eps)/(1._dp-eps*COS(z))
sint = feps*sin(sign_Ip*z)/(1._dp-eps*COS(z))
dchidr = -SIN(z)/feps**2
dchidt = sign_Ip*feps/(1._dp+eps*cost)
B = SQRT(1._dp+(eps/qbar)**2)/(1._dp+eps*cost)
dBdt = eps*sint*B/(1._dp+eps*cost)
! metric in r,chi,Phi
g11 = 1._dp
g22 = dchidr**2 + dchidt**2/eps**2
g12 = dchidr
g33 = 1._dp/(1._dp+eps*cost)**2
! magnetic field derivatives in
dBdchi=dBdt/dchidt
dBdr_c=(dBdr-g12*dBdt/dchidt)/g11
! metric in x,y,z
gxx(iz,eo) = dxdr_**2*g11
gyy(iz,eo) = (Cy*q0)**2 * (fac2*z)**2*g11 &
+ 2._dp*fac2*z*g12 &
+ Cy**2*g33
gxy(iz,eo) = dxdr_*Cy*sign_Ip*q0*(fac2*z*g11 + g12)
gxz(iz,eo) = dxdr_*g12
gyz(iz,eo) = Cy*q0*sign_Ip*(fac2*z*g12 + g22)
gzz(iz,eo) = g22
! Jacobian
Jacobian(iz,eo) = Cxy*abs(q0)*(1._dp+eps*cost)**2
! Background equilibrium magnetic field
hatB (iz,eo) = B
hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
gradxB(iz,eo) = dBdr_c/dxdr_ ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = dBdchi! Gene put a factor hatB or 1/hatR in this
! Cylindrical coordinates derivatives
dxdR(iz,eo) = cost
dxdZ(iz,eo) = sint
hatR(iz,eo) = 1._dp + eps*cost
hatZ(iz,eo) = 1._dp + eps*sint
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = X
Zc (iz,eo) = hatZ(iz,eo)
Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo)
Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo) ! Kx
Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo) ! Ky
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
! Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - (COS(z) + (shear*z - alpha_MHD*SIN(z))* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
Ckxky(iky, ikx, iz,eo) = (-sint*kx - cost*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_circular_geometry_GENE
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_zpinch_geometry
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z, kx, ky, alpha_MHD
alpha_MHD = 0._dp
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = 0._dp
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp
gyz(iz,eo) = 0._dp
gzz(iz,eo) = 1._dp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
hatR(iz,eo) = 1._dp
hatZ(iz,eo) = 1._dp
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
hatB (iz,eo) = 1._dp
hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_zpinch_geometry
!
!--------------------------------------------------------------------------------
!
subroutine eval_1D_geometry
! evaluate 1D perp geometry model
implicit none
REAL(dp) :: z, kx, ky
parity: DO eo = 0,1
zloop: DO iz = izs,ize
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp
! Relative strengh of radius
hatR(iz,eo) = 1._dp
! Jacobian
Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
hatB(iz,eo) = 1._dp
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(ikx, iky, iz,eo) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_1D_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE set_ikx_zBC_map
IMPLICIT NONE
REAL :: shift, kx_shift
! For periodic CHI BC or 0 dirichlet
LOGICAL :: PERIODIC_CHI_BC = .TRUE.
ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
! No periodic connection for extension of the domain
IF(Nexc .GT. 1) PERIODIC_CHI_BC = .TRUE.
!! No shear case (simple id mapping)
!3 | 1 2 3 4 5 6 | ky = 3 dky
!2 ky | 1 2 3 4 5 6 | ky = 2 dky
!1 A | 1 2 3 4 5 6 | ky = 1 dky
!0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
!(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=free)
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
ikx_zBC_L(iky,ikx) = ikx
ikx_zBC_R(iky,ikx) = ikx
ENDDO
ENDDO
! Modify connection map only at border of z
IF(SHEARED) THEN
! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (even NZ)
!3 | 4 x x x 2 3 | ky = 3 dky
!2 ky | 3 4 x x 1 2 | ky = 2 dky
!1 A | 2 3 4 x 6 1 | ky = 1 dky
!0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
!kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (ODD NZ)
!3 | x x x 2 3 | ky = 3 dky
!2 ky | 3 x x 1 2 | ky = 2 dky
!1 A | 2 3 x 5 1 | ky = 1 dky
!0 | -> kx | 1____2____3____4____5 | ky = 0 dky
!kx = 0 0.1 0.2 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
IF(contains_zmax) THEN ! Check if the process is at the end of the FT
DO iky = ikys,ikye
shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) + shift
! We use EPSILON() to treat perfect equality case
IF( ((kx_shift-EPSILON(kx_shift)) .GT. kx_max) .AND. (.NOT. PERIODIC_CHI_BC) )THEN ! outside of the frequ domain
ikx_zBC_R(iky,ikx) = -99
ELSE
ikx_zBC_R(iky,ikx) = ikx+(iky-1)*Nexc
IF( ikx_zBC_R(iky,ikx) .GT. Nkx ) &
ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - Nkx
ENDIF
ENDDO
ENDDO
ENDIF
! connection map BC of the LEFT boundary (z=-pi*Npol)
!3 | x 5 6 1 x x | ky = 3 dky
!2 ky | 5 6 1 2 x x | ky = 2 dky
!1 A | 6 1 2 3 x 5 | ky = 1 dky
!0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
!(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
IF(contains_zmin) THEN ! Check if the process is at the start of the FT
DO iky = ikys,ikye
shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) - shift
! We use EPSILON() to treat perfect equality case
IF( ((kx_shift+EPSILON(kx_shift)) .LT. kx_min) .AND. (.NOT. PERIODIC_CHI_BC) ) THEN ! outside of the frequ domain
ikx_zBC_L(iky,ikx) = -99
ELSE
ikx_zBC_L(iky,ikx) = ikx-(iky-1)*Nexc
IF( ikx_zBC_L(iky,ikx) .LT. 1 ) &
ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + Nkx
ENDIF
ENDDO
ENDDO
ENDIF
ELSE
ENDIF
END SUBROUTINE set_ikx_zBC_map
!
!--------------------------------------------------------------------------------
!
SUBROUTINE geometry_allocate_mem
! Curvature and geometry
CALL allocate_array( Ckxky, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
CALL allocate_array( Jacobian,izgs,izge, 0,1)
CALL allocate_array( gxx,izgs,izge, 0,1)
CALL allocate_array( gxy,izgs,izge, 0,1)
CALL allocate_array( gxz,izgs,izge, 0,1)
CALL allocate_array( gyy,izgs,izge, 0,1)
CALL allocate_array( gyz,izgs,izge, 0,1)
CALL allocate_array( gzz,izgs,izge, 0,1)
CALL allocate_array( gradxB,izgs,izge, 0,1)
CALL allocate_array( gradyB,izgs,izge, 0,1)
CALL allocate_array( gradzB,izgs,izge, 0,1)
CALL allocate_array( hatB,izgs,izge, 0,1)
CALL allocate_array( hatB_NL,izgs,izge, 0,1)
CALL allocate_array( hatR,izgs,izge, 0,1)
CALL allocate_array( hatZ,izgs,izge, 0,1)
CALL allocate_array( Rc,izgs,izge, 0,1)
CALL allocate_array( phic,izgs,izge, 0,1)
CALL allocate_array( Zc,izgs,izge, 0,1)
CALL allocate_array( dxdR,izgs,izge, 0,1)
CALL allocate_array( dxdZ,izgs,izge, 0,1)
call allocate_array(gradz_coeff,izgs,izge, 0,1)
CALL allocate_array( kparray, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
END SUBROUTINE geometry_allocate_mem
SUBROUTINE geometry_outputinputs(fidres, str)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach
USE prec_const
IMPLICIT NONE
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
CALL attach(fidres, TRIM(str),"geometry", geom)
CALL attach(fidres, TRIM(str), "q0", q0)
CALL attach(fidres, TRIM(str), "shear", shear)
CALL attach(fidres, TRIM(str), "eps", eps)
END SUBROUTINE geometry_outputinputs
end module geometry