2078172070

Committed 03 Oct 2025 09:22AM UTC coverage: 10.835% (+0.03%) from 10.806%

Build # 2078172070

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Merged in experimental/user_controlled_output_base_name (pull request #1184)

Run Details

5049 of 46599 relevant lines covered (10.83%)

119786.99 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

1.66

/src/theta_grid_file.f90

#include "unused_dummy.inc"
!> Use output from [[rungridgen]]
module theta_grid_file
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  use constants, only: run_name_size
  implicit none

  private

  public :: init_theta_grid_file, finish_theta_grid_file
  public :: check_theta_grid_file, wnml_theta_grid_file
  public :: file_get_sizes, file_get_grids
  public :: check_theta_grid_file_nc
  public :: file_nc_get_sizes, file_nc_get_grids
  public :: ntheta, nperiod, ntgrid, nbset

  public :: theta_grid_file_config_type
  public :: set_theta_grid_file_config
  public :: get_theta_grid_file_config

  character (len = run_name_size) :: gridout_file
  real :: shat_input, drhodpsi_input, kxfac_input, qval_input
  logical :: no_geo_info = .false.
  integer :: ntheta, nperiod, ntgrid, nbset
  logical :: initialized = .false.

  !> Used to represent the input configuration of theta_grid
  type, extends(abstract_config_type) :: theta_grid_file_config_type
     ! namelist : theta_grid_file_knobs
     ! indexed : false
     !> Name of file with output from [[rungridgen]].
     character(len = run_name_size) :: gridout_file = "grid.out"
     !> If false, read `Rplot`, `Rprime`, `Zplot`, `Zprime`, `aplot`, `aprime`
     !> from [[theta_grid_file_knobs:gridout_file]].
     logical :: no_geo_info = .false.
#include "theta_grid_file_overrides_and_bound_auto_gen.inc"
  end type theta_grid_file_config_type

  type(theta_grid_file_config_type) :: theta_grid_file_config

contains

  !> FIXME : Add documentation
  subroutine wnml_theta_grid_file(unit)
    implicit none
    integer, intent(in) :: unit
    call theta_grid_file_config%write(unit)
  end subroutine wnml_theta_grid_file

  !> FIXME : Add documentation
  subroutine check_theta_grid_file(report_unit)
    implicit none
    integer, intent(in) :: report_unit
    integer :: i, iunit
    real :: drhodpsi, kxfac, rmaj, shat
    character (200) :: line

    write (report_unit, *)
    write (report_unit, fmt="('Equilibrium information obtained from gridgen output file:')")
    write (report_unit, fmt="(a)") trim(gridout_file)

    open (newunit=iunit, file=gridout_file, status="old", err=100)
    read (unit=iunit, fmt="(a)") line
    read (unit=iunit, fmt=*) nbset
    read (unit=iunit, fmt="(a)") line
    do i = 1, nbset
       read (unit=iunit, fmt="(a)") line
    end do

    read (unit=iunit, fmt="(a)") line
    read (unit=iunit, fmt=*) ntgrid, nperiod, ntheta, &
         drhodpsi, rmaj, shat, kxfac

    close (unit=iunit)

    write (report_unit, *)
    write (report_unit, fmt="('Limited information available:')")
    write (report_unit, *)
    write (report_unit, fmt="('nbset =     ',i5)") nbset
    write (report_unit, fmt="('ntgrid =    ',i5)") ntgrid
    write (report_unit, fmt="('ntheta =    ',i5)") ntheta
    write (report_unit, fmt="('nperiod =   ',i2)") nperiod
    write (report_unit, fmt="('drhodpsi =  ',f8.4)") drhodpsi
    write (report_unit, fmt="('R =         ',f8.4)") Rmaj
    write (report_unit, fmt="('s_hat =     ',f8.4)") shat
    write (report_unit, fmt="('kxfac =     ',f8.4)") kxfac

    write (report_unit, *)
    write (report_unit, *) 'NOTE: Regardless of the values of ntheta and nperiod'
    write (report_unit, *) '      found in the theta_grid_parameters namelist,'
    write (report_unit, *) '      this calculation will use the values listed here:'
    write (report_unit, fmt="('ntgrid =    ',i5)") ntgrid
    write (report_unit, fmt="('ntheta =    ',i5)") ntheta
    write (report_unit, *) '      These were obtained from the gridgen output file.'
    write (report_unit, *)
100 continue
  end subroutine check_theta_grid_file

  !> FIXME : Add documentation
  subroutine check_theta_grid_file_nc(report_unit)
    use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
    use gs2_io_grids, only: nc_get_grid_sizes, nc_get_grid_scalars
    use gs2_io_grids, only: nc_get_lambda_grid_size
    implicit none
    integer, intent(in) :: report_unit
    real :: drhodpsi, kxfac, rmaj, shat, qval

    write (report_unit, *)
    write (report_unit, fmt="('Equilibrium information obtained from gridgen output file:')")
    write (report_unit, fmt="(a)") trim(gridout_file)

    call nc_grid_file_open(gridout_file, "r")
    call nc_get_grid_sizes(ntheta=ntheta, ntgrid=ntgrid, nperiod=nperiod)
    call nc_get_lambda_grid_size(nbset)
    call nc_get_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)
    call nc_grid_file_close()

    write (report_unit, *)
    write (report_unit, fmt="('Limited information available:')")
    write (report_unit, *)
    write (report_unit, fmt="('nbset =     ',i5)") nbset
    write (report_unit, fmt="('ntgrid =    ',i5)") ntgrid
    write (report_unit, fmt="('ntheta =    ',i5)") ntheta
    write (report_unit, fmt="('nperiod =   ',i2)") nperiod
    write (report_unit, fmt="('drhodpsi =  ',f8.4)") drhodpsi
    write (report_unit, fmt="('R =         ',f8.4)") Rmaj
    write (report_unit, fmt="('s_hat =     ',f8.4)") shat
    write (report_unit, fmt="('kxfac =     ',f8.4)") kxfac
    write (report_unit, fmt="('qval =     ',f8.4)") qval

    write (report_unit, *)
    write (report_unit, *) 'NOTE: Regardless of the values of ntheta and nperiod'
    write (report_unit, *) '      found in the theta_grid_parameters namelist,'
    write (report_unit, *) '      this calculation will use the values listed here:'
    write (report_unit, fmt="('ntgrid =    ',i5)") ntgrid
    write (report_unit, fmt="('ntheta =    ',i5)") ntheta
    write (report_unit, *) '      These were obtained from the gridgen output file.'
    write (report_unit, *)
  end subroutine check_theta_grid_file_nc

  !> FIXME : Add documentation
  subroutine init_theta_grid_file(theta_grid_file_config_in)
    use theta_grid_params, only: init_theta_grid_params
    implicit none
    type(theta_grid_file_config_type), intent(in), optional :: theta_grid_file_config_in

    if (initialized) return
    initialized = .true.

    call init_theta_grid_params
    call read_parameters(theta_grid_file_config_in)
  end subroutine init_theta_grid_file

  !> FIXME : Add documentation
  subroutine finish_theta_grid_file
    implicit none
    initialized = .false.
  end subroutine finish_theta_grid_file

  !> FIXME : Add documentation
  subroutine read_parameters(theta_grid_file_config_in)
    implicit none
    type(theta_grid_file_config_type), intent(in), optional :: theta_grid_file_config_in

    if (present(theta_grid_file_config_in)) theta_grid_file_config = theta_grid_file_config_in

    call theta_grid_file_config%init(name = 'theta_grid_file_knobs', requires_index = .false.)

    ! Copy out internal values into module level parameters
    associate(self => theta_grid_file_config)
#include "theta_grid_file_copy_out_auto_gen.inc"
    end associate
  end subroutine read_parameters

  !> FIXME : Add documentation
  subroutine file_get_sizes
    implicit none
    integer :: unit
    character(200) :: line
    integer :: i, ntgrid
    real :: rmaj
    open (newunit=unit, file=gridout_file, status="old")

    read (unit=unit, fmt="(a)") line
    read (unit=unit, fmt=*) nbset
    read (unit=unit, fmt="(a)") line
    do i = 1, nbset
       read (unit=unit, fmt="(a)") line
    end do

    read (unit=unit, fmt="(a)") line
    read (unit=unit, fmt=*) ntgrid, nperiod, ntheta, &
         drhodpsi_input, rmaj, shat_input, kxfac_input, qval_input

    close (unit=unit)
  end subroutine file_get_sizes

  !> FIXME : Add documentation
  subroutine file_nc_get_sizes
    use gs2_io_grids, only: nc_grid_file_open
    use gs2_io_grids, only: nc_get_grid_sizes, nc_get_lambda_grid_size
    implicit none
    call nc_grid_file_open(gridout_file, "r")
    call nc_get_grid_sizes(ntheta=ntheta, ntgrid=ntgrid, nperiod=nperiod)
    call nc_get_lambda_grid_size(nbset)
  end subroutine file_nc_get_sizes

  !> FIXME : Add documentation
  subroutine file_get_grids (nperiod, ntheta, ntgrid, nbset, theta, bset, &
       bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
       gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
       Rplot, Zplot, Rprime, Zprime, aplot, aprime, &
       shat, drhodpsi, kxfac, qval, gb_to_cv, Bpol, surfarea, dvdrhon, rhoc)
    use constants, only: pi
    use integration, only: trapezoidal_integration
    use theta_grid_params, only: rhoc_par => rhoc
    implicit none
    integer, intent (in) :: nperiod
    integer, intent (in out) :: ntheta, ntgrid, nbset
    real, dimension (-ntgrid:ntgrid), intent (out) :: theta
    real, dimension (nbset), intent (out) :: bset
    real, dimension (-ntgrid:ntgrid), intent (out) :: &
         bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
         gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
         Rplot, Zplot, Rprime, Zprime, aplot, aprime, Bpol
    real, intent (out) :: shat, drhodpsi, kxfac, qval, surfarea, dvdrhon, rhoc
    logical, intent (in) :: gb_to_cv
    integer :: unit
    character(200) :: line
    integer :: i
    logical, save :: first = .true.
    !Not used because we don't call gridgen_get_grids for file
    UNUSED_DUMMY(nperiod); UNUSED_DUMMY(ntheta)

    !<DD> Should jacob also be provided by this routine?

    !<DD> NOTE: Not currently settin Bpol here. This is used in flux calculations.
    !If not set then results will be funny. Add a warning message and set to zero
    !for now.
    if (first) then
       write(*,'("WARNING: When using file_get_grids, Bpol does not have a correct definition --> Currently just setting to zero, may impact some diagnostics")')
       Bpol = 0.
       first = .false.
    end if

    shat = shat_input
    drhodpsi = drhodpsi_input
    kxfac = kxfac_input
    qval = qval_input

    open (newunit=unit, file=gridout_file, status="old")
    read (unit=unit, fmt="(a)") line
    read (unit=unit, fmt="(a)") line
    read (unit=unit, fmt="(a)") line
    do i = 1, nbset
       read (unit=unit, fmt=*) bset(i) ! actually alambda
    end do
    bset = 1.0/bset ! switch alambda to bset

    read (unit=unit, fmt="(a)") line
    read (unit=unit, fmt="(a)") line

    read (unit=unit, fmt="(a)") line
    do i = -ntgrid, ntgrid
       read (unit=unit, fmt=*) gbdrift(i), gradpar(i), grho(i)
    end do

    read (unit=unit, fmt="(a)") line
    do i = -ntgrid, ntgrid
       read (unit=unit, fmt=*) cvdrift(i), gds2(i), bmag(i), theta(i)
    end do

    read (unit=unit, fmt="(a)") line
    do i = -ntgrid, ntgrid
       read (unit=unit, fmt=*) gds21(i), gds22(i)
    end do

    ! TMP UNTIL WORK OUT HOW TO GET FROM FILE
    gds23 = 0. ; gds24 = 0. ; gds24_noq = 0.

    ! TMP UNTIL FIGURE OUT HOW TO WORK WITH FILE -- MAB
    ! set coriolis drift to zero
    cdrift = 0. ; cdrift0 = 0.

    read (unit=unit, fmt="(a)") line
    do i = -ntgrid, ntgrid
       read (unit=unit, fmt=*) cvdrift0(i), gbdrift0(i)
    end do

    if (gb_to_cv) then
       do i =-ntgrid,ntgrid
          gbdrift(i) = cvdrift(i)
          gbdrift0(i) = cvdrift0(i)
       end do
    end if

    ! Possibly crude approximations for file
    ! Note jacob is usually defined as jacob = 1.0/(drhodpsi*gradpar*bmag)
    ! but isn't yet stored here.
    surfarea = 2 * pi * trapezoidal_integration(theta, grho / (drhodpsi*gradpar*bmag))
    dvdrhon = 2 * pi * trapezoidal_integration(theta, 1.0 / (drhodpsi*gradpar*bmag))

    ! As the eqfile doesn't specify rhoc we simply set from the value
    ! from theta_grid_params here.
    rhoc = rhoc_par

    if (.not. no_geo_info) then

       read (unit=unit, fmt="(a)",err=100) line
       do i = -ntgrid, ntgrid
          read (unit=unit, fmt=*, err=100) Rplot(i), Rprime(i)
       end do

       read (unit=unit, fmt="(a)",err=100) line
       do i = -ntgrid, ntgrid
          read (unit=unit, fmt=*, err=100) Zplot(i), Zprime(i)
       end do

       read (unit=unit, fmt="(a)",err=100) line
       do i = -ntgrid, ntgrid
          read (unit=unit, fmt=*, err=100) aplot(i), aprime(i)
       end do

    end if

    close (unit=unit)

    return

100 continue
    write(6,*) 'Error reading Rplot etc. setting to dummy values.'
! dummy values for backward compatibility
    Rplot = 1. ; Rprime = 0.
    Zplot = 1. ; Zprime = 0.
    aplot = 1. ; aprime = 0.

    close (unit=unit)

  end subroutine file_get_grids

  !> FIXME : Add documentation
  subroutine file_nc_get_grids (nperiod, ntheta, ntgrid, nbset, theta, bset, &
       bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
       gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
       Rplot, Zplot, Rprime, Zprime, aplot, aprime, &
       shat, drhodpsi, kxfac, qval, gb_to_cv, Bpol, surfarea, dvdrhon, rhoc)
    use constants, only: pi
    use integration, only: trapezoidal_integration
    use theta_grid_params, only: rhoc_par => rhoc
    use gs2_io_grids, only: nc_grid_file_close
    use gs2_io_grids, only: nc_get_grids, nc_get_grid_scalars
    use gs2_io_grids, only: nc_get_lambda_grid
    implicit none
    integer, intent (in) :: nperiod
    integer, intent (in out) :: ntheta, ntgrid, nbset
    real, dimension (-ntgrid:ntgrid), intent (out) :: theta
    real, dimension (nbset), intent (out) :: bset
    real, dimension (-ntgrid:ntgrid), intent (out) :: &
         bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
         gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
         Rplot, Zplot, Rprime, Zprime, aplot, aprime, Bpol
    real, intent (out) :: shat, drhodpsi, kxfac, qval, surfarea, dvdrhon, rhoc
    real :: rmaj
    logical, intent (in) :: gb_to_cv
    !Not used because we don't call gridgen_get_grids for file
    UNUSED_DUMMY(nperiod); UNUSED_DUMMY(ntheta)
    Rplot = 1.; Rprime = 0.
    Zplot = 1.; Zprime = 0.
    aplot = 1.; aprime = 0.

    call nc_get_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)

    call nc_get_grids(ntgrid, &
       bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, &
       gds2, gds21, gds22, grho, theta, &
       cdrift=cdrift, cdrift0=cdrift0, &
       Rplot=Rplot, Rprime=Rprime, Zplot=Zplot, Zprime=Zprime, aplot=aplot, aprime=aprime, &
       Bpol=Bpol)

    call nc_get_lambda_grid(nbset, bset)

    call nc_grid_file_close()

    ! switch alambda to bset
    bset = 1.0/bset

    ! TMP UNTIL WORK OUT HOW TO GET FROM FILE
    gds23 = 0. ; gds24 = 0. ; gds24_noq = 0.

    ! TMP UNTIL FIGURE OUT HOW TO WORK WITH FILE -- MAB
    ! set coriolis drift to zero
    cdrift = 0. ; cdrift0 = 0.

    if (gb_to_cv) then
       gbdrift = cvdrift
       gbdrift0 = cvdrift0
    end if

    ! Possibly crude approximations for file
    ! Note jacob is usually defined as jacob = 1.0/(drhodpsi*gradpar*bmag)
    ! but isn't yet stored here.
    surfarea = 2 * pi * trapezoidal_integration(theta, grho / (drhodpsi*gradpar*bmag))
    dvdrhon = 2 * pi * trapezoidal_integration(theta, 1.0 / (drhodpsi*gradpar*bmag))

    ! As the eqfile doesn't specify rhoc we simply set from the value
    ! from theta_grid_params here.
    rhoc = rhoc_par

  end subroutine file_nc_get_grids

#include "theta_grid_file_auto_gen.inc"
end module theta_grid_file

1	#include "unused_dummy.inc"
2	!> Use output from [[rungridgen]]
3	module theta_grid_file
4	use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
5	use constants, only: run_name_size
6	implicit none
7
8	private
9
10	public :: init_theta_grid_file, finish_theta_grid_file
11	public :: check_theta_grid_file, wnml_theta_grid_file
12	public :: file_get_sizes, file_get_grids
13	public :: check_theta_grid_file_nc
14	public :: file_nc_get_sizes, file_nc_get_grids
15	public :: ntheta, nperiod, ntgrid, nbset
16
17	public :: theta_grid_file_config_type
18	public :: set_theta_grid_file_config
19	public :: get_theta_grid_file_config
20
21	character (len = run_name_size) :: gridout_file
22	real :: shat_input, drhodpsi_input, kxfac_input, qval_input
23	logical :: no_geo_info = .false.
24	integer :: ntheta, nperiod, ntgrid, nbset
25	logical :: initialized = .false.
26
27	!> Used to represent the input configuration of theta_grid
28	type, extends(abstract_config_type) :: theta_grid_file_config_type
29	! namelist : theta_grid_file_knobs
30	! indexed : false
31	!> Name of file with output from [[rungridgen]].
32	character(len = run_name_size) :: gridout_file = "grid.out"
33	!> If false, read `Rplot`, `Rprime`, `Zplot`, `Zprime`, `aplot`, `aprime`
34	!> from [[theta_grid_file_knobs:gridout_file]].
35	logical :: no_geo_info = .false.
36	#include "theta_grid_file_overrides_and_bound_auto_gen.inc"
37	end type theta_grid_file_config_type
38
39	type(theta_grid_file_config_type) :: theta_grid_file_config
40
41	contains
42
43	!> FIXME : Add documentation
44	subroutine wnml_theta_grid_file(unit)	×
45	implicit none
46	integer, intent(in) :: unit
47	call theta_grid_file_config%write(unit)	×
48	end subroutine wnml_theta_grid_file	×
49
50	!> FIXME : Add documentation
51	subroutine check_theta_grid_file(report_unit)	×
52	implicit none
53	integer, intent(in) :: report_unit
54	integer :: i, iunit
55	real :: drhodpsi, kxfac, rmaj, shat
56	character (200) :: line
57
58	write (report_unit, *)	×
59	write (report_unit, fmt="('Equilibrium information obtained from gridgen output file:')")	×
60	write (report_unit, fmt="(a)") trim(gridout_file)	×
61
62	open (newunit=iunit, file=gridout_file, status="old", err=100)	×
63	read (unit=iunit, fmt="(a)") line	×
64	read (unit=iunit, fmt=*) nbset	×
65	read (unit=iunit, fmt="(a)") line	×
66	do i = 1, nbset	×
67	read (unit=iunit, fmt="(a)") line	×
68	end do
69
70	read (unit=iunit, fmt="(a)") line	×
71	read (unit=iunit, fmt=*) ntgrid, nperiod, ntheta, &	×
72	drhodpsi, rmaj, shat, kxfac	×
73
74	close (unit=iunit)	×
75
76	write (report_unit, *)	×
77	write (report_unit, fmt="('Limited information available:')")	×
78	write (report_unit, *)	×
79	write (report_unit, fmt="('nbset = ',i5)") nbset	×
80	write (report_unit, fmt="('ntgrid = ',i5)") ntgrid	×
81	write (report_unit, fmt="('ntheta = ',i5)") ntheta	×
82	write (report_unit, fmt="('nperiod = ',i2)") nperiod	×
83	write (report_unit, fmt="('drhodpsi = ',f8.4)") drhodpsi	×
84	write (report_unit, fmt="('R = ',f8.4)") Rmaj	×
85	write (report_unit, fmt="('s_hat = ',f8.4)") shat	×
86	write (report_unit, fmt="('kxfac = ',f8.4)") kxfac	×
87
88	write (report_unit, *)	×
89	write (report_unit, *) 'NOTE: Regardless of the values of ntheta and nperiod'	×
90	write (report_unit, *) ' found in the theta_grid_parameters namelist,'	×
91	write (report_unit, *) ' this calculation will use the values listed here:'	×
92	write (report_unit, fmt="('ntgrid = ',i5)") ntgrid	×
93	write (report_unit, fmt="('ntheta = ',i5)") ntheta	×
94	write (report_unit, *) ' These were obtained from the gridgen output file.'	×
95	write (report_unit, *)	×
96	100 continue
97	end subroutine check_theta_grid_file	×
98
99	!> FIXME : Add documentation
100	subroutine check_theta_grid_file_nc(report_unit)	×
101	use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
102	use gs2_io_grids, only: nc_get_grid_sizes, nc_get_grid_scalars
103	use gs2_io_grids, only: nc_get_lambda_grid_size
104	implicit none
105	integer, intent(in) :: report_unit
106	real :: drhodpsi, kxfac, rmaj, shat, qval
107
108	write (report_unit, *)	×
109	write (report_unit, fmt="('Equilibrium information obtained from gridgen output file:')")	×
110	write (report_unit, fmt="(a)") trim(gridout_file)	×
111
112	call nc_grid_file_open(gridout_file, "r")	×
113	call nc_get_grid_sizes(ntheta=ntheta, ntgrid=ntgrid, nperiod=nperiod)	×
114	call nc_get_lambda_grid_size(nbset)	×
115	call nc_get_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)	×
116	call nc_grid_file_close()	×
117
118	write (report_unit, *)	×
119	write (report_unit, fmt="('Limited information available:')")	×
120	write (report_unit, *)	×
121	write (report_unit, fmt="('nbset = ',i5)") nbset	×
122	write (report_unit, fmt="('ntgrid = ',i5)") ntgrid	×
123	write (report_unit, fmt="('ntheta = ',i5)") ntheta	×
124	write (report_unit, fmt="('nperiod = ',i2)") nperiod	×
125	write (report_unit, fmt="('drhodpsi = ',f8.4)") drhodpsi	×
126	write (report_unit, fmt="('R = ',f8.4)") Rmaj	×
127	write (report_unit, fmt="('s_hat = ',f8.4)") shat	×
128	write (report_unit, fmt="('kxfac = ',f8.4)") kxfac	×
129	write (report_unit, fmt="('qval = ',f8.4)") qval	×
130
131	write (report_unit, *)	×
132	write (report_unit, *) 'NOTE: Regardless of the values of ntheta and nperiod'	×
133	write (report_unit, *) ' found in the theta_grid_parameters namelist,'	×
134	write (report_unit, *) ' this calculation will use the values listed here:'	×
135	write (report_unit, fmt="('ntgrid = ',i5)") ntgrid	×
136	write (report_unit, fmt="('ntheta = ',i5)") ntheta	×
137	write (report_unit, *) ' These were obtained from the gridgen output file.'	×
138	write (report_unit, *)	×
139	end subroutine check_theta_grid_file_nc	×
140
141	!> FIXME : Add documentation
142	subroutine init_theta_grid_file(theta_grid_file_config_in)	×
143	use theta_grid_params, only: init_theta_grid_params
144	implicit none
145	type(theta_grid_file_config_type), intent(in), optional :: theta_grid_file_config_in
146
147	if (initialized) return	×
148	initialized = .true.	×
149
150	call init_theta_grid_params	×
151	call read_parameters(theta_grid_file_config_in)	×
152	end subroutine init_theta_grid_file
153
154	!> FIXME : Add documentation
155	subroutine finish_theta_grid_file	28✔
156	implicit none
157	initialized = .false.	28✔
158	end subroutine finish_theta_grid_file	28✔
159
160	!> FIXME : Add documentation
161	subroutine read_parameters(theta_grid_file_config_in)	×
162	implicit none
163	type(theta_grid_file_config_type), intent(in), optional :: theta_grid_file_config_in
164
165	if (present(theta_grid_file_config_in)) theta_grid_file_config = theta_grid_file_config_in	×
166
167	call theta_grid_file_config%init(name = 'theta_grid_file_knobs', requires_index = .false.)	×
168
169	! Copy out internal values into module level parameters
170	associate(self => theta_grid_file_config)
171	#include "theta_grid_file_copy_out_auto_gen.inc"
172	end associate
173	end subroutine read_parameters	×
174
175	!> FIXME : Add documentation
176	subroutine file_get_sizes	×
177	implicit none
178	integer :: unit
179	character(200) :: line
180	integer :: i, ntgrid
181	real :: rmaj
182	open (newunit=unit, file=gridout_file, status="old")	×
183
184	read (unit=unit, fmt="(a)") line	×
185	read (unit=unit, fmt=*) nbset	×
186	read (unit=unit, fmt="(a)") line	×
187	do i = 1, nbset	×
188	read (unit=unit, fmt="(a)") line	×
189	end do
190
191	read (unit=unit, fmt="(a)") line	×
192	read (unit=unit, fmt=*) ntgrid, nperiod, ntheta, &	×
193	drhodpsi_input, rmaj, shat_input, kxfac_input, qval_input	×
194
195	close (unit=unit)	×
196	end subroutine file_get_sizes	×
197
198	!> FIXME : Add documentation
199	subroutine file_nc_get_sizes	×
200	use gs2_io_grids, only: nc_grid_file_open
201	use gs2_io_grids, only: nc_get_grid_sizes, nc_get_lambda_grid_size
202	implicit none
203	call nc_grid_file_open(gridout_file, "r")	×
204	call nc_get_grid_sizes(ntheta=ntheta, ntgrid=ntgrid, nperiod=nperiod)	×
205	call nc_get_lambda_grid_size(nbset)	×
206	end subroutine file_nc_get_sizes	×
207
208	!> FIXME : Add documentation
209	subroutine file_get_grids (nperiod, ntheta, ntgrid, nbset, theta, bset, &	×
210	bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &	×
211	gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &	×
212	Rplot, Zplot, Rprime, Zprime, aplot, aprime, &	×
213	shat, drhodpsi, kxfac, qval, gb_to_cv, Bpol, surfarea, dvdrhon, rhoc)	×
214	use constants, only: pi
215	use integration, only: trapezoidal_integration
216	use theta_grid_params, only: rhoc_par => rhoc
217	implicit none
218	integer, intent (in) :: nperiod
219	integer, intent (in out) :: ntheta, ntgrid, nbset
220	real, dimension (-ntgrid:ntgrid), intent (out) :: theta
221	real, dimension (nbset), intent (out) :: bset
222	real, dimension (-ntgrid:ntgrid), intent (out) :: &
223	bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
224	gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
225	Rplot, Zplot, Rprime, Zprime, aplot, aprime, Bpol
226	real, intent (out) :: shat, drhodpsi, kxfac, qval, surfarea, dvdrhon, rhoc
227	logical, intent (in) :: gb_to_cv
228	integer :: unit
229	character(200) :: line
230	integer :: i
231	logical, save :: first = .true.
232	!Not used because we don't call gridgen_get_grids for file
233	UNUSED_DUMMY(nperiod); UNUSED_DUMMY(ntheta)	×
234
235	!<DD> Should jacob also be provided by this routine?
236
237	!<DD> NOTE: Not currently settin Bpol here. This is used in flux calculations.
238	!If not set then results will be funny. Add a warning message and set to zero
239	!for now.
240	if (first) then	×
241	write(*,'("WARNING: When using file_get_grids, Bpol does not have a correct definition --> Currently just setting to zero, may impact some diagnostics")')	×
242	Bpol = 0.	×
243	first = .false.	×
244	end if
245
246	shat = shat_input	×
247	drhodpsi = drhodpsi_input	×
248	kxfac = kxfac_input	×
249	qval = qval_input	×
250
251	open (newunit=unit, file=gridout_file, status="old")	×
252	read (unit=unit, fmt="(a)") line	×
253	read (unit=unit, fmt="(a)") line	×
254	read (unit=unit, fmt="(a)") line	×
255	do i = 1, nbset	×
256	read (unit=unit, fmt=*) bset(i) ! actually alambda	×
257	end do
258	bset = 1.0/bset ! switch alambda to bset	×
259
260	read (unit=unit, fmt="(a)") line	×
261	read (unit=unit, fmt="(a)") line	×
262
263	read (unit=unit, fmt="(a)") line	×
264	do i = -ntgrid, ntgrid	×
265	read (unit=unit, fmt=*) gbdrift(i), gradpar(i), grho(i)	×
266	end do
267
268	read (unit=unit, fmt="(a)") line	×
269	do i = -ntgrid, ntgrid	×
270	read (unit=unit, fmt=*) cvdrift(i), gds2(i), bmag(i), theta(i)	×
271	end do
272
273	read (unit=unit, fmt="(a)") line	×
274	do i = -ntgrid, ntgrid	×
275	read (unit=unit, fmt=*) gds21(i), gds22(i)	×
276	end do
277
278	! TMP UNTIL WORK OUT HOW TO GET FROM FILE
279	gds23 = 0. ; gds24 = 0. ; gds24_noq = 0.	×
280
281	! TMP UNTIL FIGURE OUT HOW TO WORK WITH FILE -- MAB
282	! set coriolis drift to zero
283	cdrift = 0. ; cdrift0 = 0.	×
284
285	read (unit=unit, fmt="(a)") line	×
286	do i = -ntgrid, ntgrid	×
287	read (unit=unit, fmt=*) cvdrift0(i), gbdrift0(i)	×
288	end do
289
290	if (gb_to_cv) then	×
291	do i =-ntgrid,ntgrid	×
292	gbdrift(i) = cvdrift(i)	×
293	gbdrift0(i) = cvdrift0(i)	×
294	end do
295	end if
296
297	! Possibly crude approximations for file
298	! Note jacob is usually defined as jacob = 1.0/(drhodpsigradparbmag)
299	! but isn't yet stored here.
300	surfarea = 2 * pi * trapezoidal_integration(theta, grho / (drhodpsigradparbmag))	×
301	dvdrhon = 2 * pi * trapezoidal_integration(theta, 1.0 / (drhodpsigradparbmag))	×
302
303	! As the eqfile doesn't specify rhoc we simply set from the value
304	! from theta_grid_params here.
305	rhoc = rhoc_par	×
306
307	if (.not. no_geo_info) then	×
308
309	read (unit=unit, fmt="(a)",err=100) line	×
310	do i = -ntgrid, ntgrid	×
311	read (unit=unit, fmt=*, err=100) Rplot(i), Rprime(i)	×
312	end do
313
314	read (unit=unit, fmt="(a)",err=100) line	×
315	do i = -ntgrid, ntgrid	×
316	read (unit=unit, fmt=*, err=100) Zplot(i), Zprime(i)	×
317	end do
318
319	read (unit=unit, fmt="(a)",err=100) line	×
320	do i = -ntgrid, ntgrid	×
321	read (unit=unit, fmt=*, err=100) aplot(i), aprime(i)	×
322	end do
323
324	end if
325
326	close (unit=unit)	×
327
328	return	×
329
330	100 continue
331	write(6,*) 'Error reading Rplot etc. setting to dummy values.'	×
332	! dummy values for backward compatibility
333	Rplot = 1. ; Rprime = 0.	×
334	Zplot = 1. ; Zprime = 0.	×
335	aplot = 1. ; aprime = 0.	×
336
337	close (unit=unit)	×
338
339	end subroutine file_get_grids
340
341	!> FIXME : Add documentation
342	subroutine file_nc_get_grids (nperiod, ntheta, ntgrid, nbset, theta, bset, &	×
343	bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &	×
344	gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &	×
345	Rplot, Zplot, Rprime, Zprime, aplot, aprime, &	×
346	shat, drhodpsi, kxfac, qval, gb_to_cv, Bpol, surfarea, dvdrhon, rhoc)	×
347	use constants, only: pi
348	use integration, only: trapezoidal_integration
349	use theta_grid_params, only: rhoc_par => rhoc
350	use gs2_io_grids, only: nc_grid_file_close
351	use gs2_io_grids, only: nc_get_grids, nc_get_grid_scalars
352	use gs2_io_grids, only: nc_get_lambda_grid
353	implicit none
354	integer, intent (in) :: nperiod
355	integer, intent (in out) :: ntheta, ntgrid, nbset
356	real, dimension (-ntgrid:ntgrid), intent (out) :: theta
357	real, dimension (nbset), intent (out) :: bset
358	real, dimension (-ntgrid:ntgrid), intent (out) :: &
359	bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0, &
360	gds2, gds21, gds22, gds23, gds24, gds24_noq, grho, &
361	Rplot, Zplot, Rprime, Zprime, aplot, aprime, Bpol
362	real, intent (out) :: shat, drhodpsi, kxfac, qval, surfarea, dvdrhon, rhoc
363	real :: rmaj
364	logical, intent (in) :: gb_to_cv
365	!Not used because we don't call gridgen_get_grids for file
366	UNUSED_DUMMY(nperiod); UNUSED_DUMMY(ntheta)	×
367	Rplot = 1.; Rprime = 0.	×
368	Zplot = 1.; Zprime = 0.	×
369	aplot = 1.; aprime = 0.	×
370
371	call nc_get_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)	×
372
373	call nc_get_grids(ntgrid, &
374	bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, &
375	gds2, gds21, gds22, grho, theta, &
376	cdrift=cdrift, cdrift0=cdrift0, &
377	Rplot=Rplot, Rprime=Rprime, Zplot=Zplot, Zprime=Zprime, aplot=aplot, aprime=aprime, &
378	Bpol=Bpol)	×
379
380	call nc_get_lambda_grid(nbset, bset)	×
381
382	call nc_grid_file_close()	×
383
384	! switch alambda to bset
385	bset = 1.0/bset	×
386
387	! TMP UNTIL WORK OUT HOW TO GET FROM FILE
388	gds23 = 0. ; gds24 = 0. ; gds24_noq = 0.	×
389
390	! TMP UNTIL FIGURE OUT HOW TO WORK WITH FILE -- MAB
391	! set coriolis drift to zero
392	cdrift = 0. ; cdrift0 = 0.	×
393
394	if (gb_to_cv) then	×
395	gbdrift = cvdrift	×
396	gbdrift0 = cvdrift0	×
397	end if
398
399	! Possibly crude approximations for file
400	! Note jacob is usually defined as jacob = 1.0/(drhodpsigradparbmag)
401	! but isn't yet stored here.
402	surfarea = 2 * pi * trapezoidal_integration(theta, grho / (drhodpsigradparbmag))	×
403	dvdrhon = 2 * pi * trapezoidal_integration(theta, 1.0 / (drhodpsigradparbmag))	×
404
405	! As the eqfile doesn't specify rhoc we simply set from the value
406	! from theta_grid_params here.
407	rhoc = rhoc_par	×
408
409	end subroutine file_nc_get_grids	×
410
411	#include "theta_grid_file_auto_gen.inc"
412	end module theta_grid_file	×

gyrokinetics / gs2 / 2078172070

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous