istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

convert remaining subroutines to derived types

Browse Source 
- converts analytical profiles/equilibrium subroutines to derived types
- use less undecipherable and consistent names

The subroutine names have changed as follows:

    set_prfspl → set_profiles_spline
    set_prfan  → set_profiles_an
  unset_prfspl → unset_profiles_spline
  unset_prfan  → unset_profiles_an
  set_equian   → set_equil_an
  set_eqspl    → set_equil_spline
  unset_equian → unset_equil_an
  unset_eqspl  → unset_equil_spline
  unset_rhospl → unset_rho_spline
This commit is contained in:
Michele Guerini Rocco 2022-05-11 00:30:34 +02:00
parent 9fcf8e51e8
commit 63e2bf0b04
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
4 changed files with 175 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
src/coreprofiles.f90
View File

@ -1,5 +1,5 @@
module coreprofiles
use const_and_precisions, only : wp_,zero,one
use const_and_precisions, only : wp_, zero, one
implicit none
integer, save :: npp, nsfd
@ -185,25 +185,33 @@ contains
end subroutine read_profiles
subroutine read_profiles_an(filenm,te,ne,zeff,unit)
use utils, only : get_free_unit
use logger, only : log_error
subroutine read_profiles_an(filenm, data, unit)
! Reads the plasma profiles `data` in the analytical format
! from params%filenm.
! If given, the file is opened in the `unit` number.
!
! TODO: add format description
use gray_params, only : profiles_data
use utils, only : get_free_unit
use logger, only : log_error
implicit none
! arguments
character(len=*), intent(in) :: filenm
real(wp_), dimension(:), allocatable, intent(out) :: te,ne,zeff
integer, optional, intent(in) :: unit
! local variables
! subroutine arguments
character(len=*), intent(in) :: filenm
type(profiles_data), intent(out) :: data
integer, optional, intent(in) :: unit
! local variables
integer :: u
integer :: err
u = get_free_unit(unit)
if(allocated(te)) deallocate(te)
if(allocated(ne)) deallocate(ne)
if(allocated(zeff)) deallocate(zeff)
allocate(te(4),ne(3),zeff(1))
if (allocated(data%terad)) deallocate(data%terad)
if (allocated(data%derad)) deallocate(data%derad)
if (allocated(data%zfc)) deallocate(data%zfc)
allocate(data%terad(4), data%derad(3), data%zfc(1))
open(file=filenm, status='old', action='read', unit=u, iostat=err)
if (err /= 0) then
@ -212,9 +220,9 @@ contains
call exit(1)
end if
read(u,*) ne(1:3) ! dens0,aln1,aln2
read(u,*) te(1:4) ! te0,dte0,alt1,alt2
read(u,*) zeff(1) ! zeffan
read (u,*) data%derad(1:3) ! dens0, aln1, aln2
read (u,*) data%terad(1:4) ! te0, dte0, alt1, alt2
read (u,*) data%zfc(1) ! zeffan
close(u)
end subroutine read_profiles_an
@ -264,7 +272,7 @@ contains
end subroutine scale_profiles
subroutine set_prfspl(params, data)
subroutine set_profiles_spline(params, data)
! Computes splines for the plasma profiles data and stores them
! in their respective global variables, see the top of this file.
use simplespline, only : difcs
@ -361,37 +369,43 @@ contains
psdbnd=min(psdbnd,xxp)
end if
write (msg, '(a,g0.3)') 'density boundary: ψ=', psdbnd
call log_info(msg, mod="coreprofiles", proc="set_prfspl")
call log_info(msg, mod="coreprofiles", proc="set_profiles_spline")
end if
deallocate(iwrkf,wrkf,wf)
end subroutine set_prfspl
end subroutine set_profiles_spline
subroutine unset_prfspl
subroutine unset_profiles_spline
implicit none
if(allocated(psrad)) deallocate(psrad)
if(allocated(ct)) deallocate(ct)
if(allocated(cz)) deallocate(cz)
if(allocated(tfn)) deallocate(tfn)
if(allocated(cfn)) deallocate(cfn)
end subroutine unset_prfspl
if (allocated(psrad)) deallocate(psrad)
if (allocated(ct)) deallocate(ct)
if (allocated(cz)) deallocate(cz)
if (allocated(tfn)) deallocate(tfn)
if (allocated(cfn)) deallocate(cfn)
end subroutine unset_profiles_spline
subroutine set_profiles_an(data)
! Stores the analytical profiles data in their respective
! global variables, see the top of this file.
use gray_params, only : profiles_data
subroutine set_prfan(te,ne,zeff)
implicit none
REAL(wp_), dimension(:), intent(in) :: te,ne,zeff
te0=te(1)
dte0=te(2)
alt1=te(3)
alt2=te(4)
dens0=ne(1)
aln1=ne(2)
aln2=ne(3)
zeffan=zeff(1)
! subroutine arguments
type(profiles_data), intent(in) :: data
psdbnd=1.0_wp_
end subroutine set_prfan
te0 = data%terad(1)
dte0 = data%terad(2)
alt1 = data%terad(3)
alt2 = data%terad(4)
dens0 = data%derad(1)
aln1 = data%derad(2)
aln2 = data%derad(3)
zeffan = data%zfc(1)
psdbnd = one
end subroutine set_profiles_an
end module coreprofiles

173
src/equilibrium.f90
View File

@ -160,20 +160,28 @@ contains
end subroutine read_eqdsk
subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit)
use utils, only : get_free_unit
use logger, only : log_error
subroutine read_equil_an(filenm, ipass, data, unit)
! Reads the MHD equilibrium `data` in the analytical format
! from params%filenm.
! If given, the file is opened in the `unit` number.
!
! TODO: add format description
use gray_params, only : equilibrium_data
use utils, only : get_free_unit
use logger, only : log_error
implicit none
! arguments
character(len=*), intent(in) :: filenm
integer, intent(in) :: ipass
integer, optional, intent(in) :: unit
real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim
! local variables
! subroutine arguments
character(len=*), intent(in) :: filenm
integer, intent(in) :: ipass
type(equilibrium_data), intent(out) :: data
integer, optional, intent(in) :: unit
! local variables
integer :: i, u, nlim
integer :: err
real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen
real(wp_) :: rr0m, zr0m, rpam, b0, q0, qa, alq
u = get_free_unit(unit)
@ -184,48 +192,37 @@ contains
call exit(1)
end if
read(u,*) rr0m,zr0m,rpam
read(u,*) b0
read(u,*) q0,qa,alq
! rcen=rr0m
! btrcen=b0
! b0=isgnbphi*b0*factb
! rvac=rr0m
! rax=rr0m
! zax=zr0m
! zbmin=(zr0-rpa)/1.0e2_wp_
! zbmax=(zr0+rpa)/1.0e2_wp_
if(allocated(rv)) deallocate(rv)
if(allocated(zv)) deallocate(zv)
if(allocated(fpol)) deallocate(fpol)
if(allocated(q)) deallocate(q)
allocate(rv(2),zv(1),fpol(1),q(3))
rv(1)=rr0m
rv(2)=rpam
zv(1)=zr0m
fpol(1)=b0*rr0m
q(1)=q0
q(2)=qa
q(3)=alq
read(u, *) rr0m, zr0m, rpam
read(u, *) b0
read(u, *) q0, qa, alq
if(ipass.ge.2) then
if(allocated(data%rv)) deallocate(data%rv)
if(allocated(data%zv)) deallocate(data%zv)
if(allocated(data%fpol)) deallocate(data%fpol)
if(allocated(data%qpsi)) deallocate(data%qpsi)
allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3))
! get size of boundary and limiter arrays and allocate them
data%rv = [rr0m, rpam]
data%zv = [zr0m]
data%fpol = [b0*rr0m]
data%qpsi = [q0, qa, alq]
if(ipass >= 2) then
! get size of boundary and limiter arrays and allocate them
read (u,*) nlim
if (allocated(rlim)) deallocate(rlim)
if (allocated(zlim)) deallocate(zlim)
! store boundary and limiter data
if(nlim>0) then
allocate(rlim(nlim),zlim(nlim))
read(u,*) (rlim(i),zlim(i),i=1,nlim)
end if
if (allocated(data%rlim)) deallocate(data%rlim)
if (allocated(data%zlim)) deallocate(data%zlim)
end if
! store boundary and limiter data
if(nlim > 0) then
allocate(data%rlim(nlim), data%zlim(nlim))
read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
end if
end if
close(u)
end subroutine read_equil_an
subroutine change_cocos(data, cocosin, cocosout, error)
! Convert the MHD equilibrium data from one coordinate convention
@ -345,7 +342,7 @@ contains
end subroutine eq_scal
subroutine set_eqspl(params, data)
subroutine set_equil_spline(params, data)
! Computes splines for the MHD equilibrium data and stores them
! in their respective global variables, see the top of this file.
use const_and_precisions, only : zero, one
@ -547,7 +544,7 @@ contains
write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
call log_info(msg, mod='equilibrium', proc='set_eqspl')
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
! search for X-point if params%ixp /= 0
@ -558,7 +555,7 @@ contains
if(psinxptmp/=-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
call log_info(msg, mod='equilibrium', proc='set_eqspl')
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
zbinf=z1
psinop=psinoptmp
@ -573,7 +570,7 @@ contains
if(psinxptmp.ne.-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
'r', r1, 'z', z1, 'ψ', psinxptmp
call log_info(msg, mod='equilibrium', proc='set_eqspl')
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
zbsup=z1
psinop=psinoptmp
@ -601,7 +598,7 @@ contains
rbinf=r1
write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
'r', rbinf, rbsup, 'z', zbinf, zbsup
call log_info(msg, mod='equilibrium', proc='set_eqspl')
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
end if
! Save Bt value on axis (required in flux_average and used in Jcd def)
@ -613,13 +610,13 @@ contains
rcen = data%rvac
write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis
call log_info(msg, mod='equilibrium', proc='set_eqspl')
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
! Compute rho_pol/rho_tor mapping based on input q profile
call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn)
call set_rhospl(sqrt(data%psinr),rhotn)
end subroutine set_eqspl
end subroutine set_equil_spline
subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
@ -706,19 +703,33 @@ contains
end subroutine setqphi_num
subroutine set_equil_an(data, n)
! Computes the analytical equilibrium data and stores them
! in their respective global variables, see the top of this file.
use const_and_precisions, only : pi, zero, one
use gray_params, only : equilibrium_data
subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n)
use const_and_precisions, only : pi,zero,one
implicit none
! arguments
real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp
integer, intent(in), optional :: n
! local variables
! subroutine arguments
type(equilibrium_data), intent(in) :: data
integer, optional, intent(in) :: n
! local variables
integer, parameter :: nqdef=101
integer :: i
real(wp_) :: dr,fq0,fq1,qq,res,rn
real(wp_) :: dr, fq0, fq1, qq, res, rn
real(wp_) :: rax, zax, a, bax, qax, q1, qexp
real(wp_), dimension(:), allocatable :: rhotn,rhopn
rax = data%rv(1)
zax = data%zv(1)
a = data%rv(2)
bax = data%fpol(1) / rax
qax = data%qpsi(1)
q1 = data%qpsi(2)
qexp = data%qpsi(3)
btaxis=bax
rmaxis=rax
zmaxis=zax
@ -767,8 +778,7 @@ contains
rhopn=sqrt(psinr)
call set_rhospl(rhopn,rhotn)
end subroutine set_equian
end subroutine set_equil_an
subroutine set_rhospl(rhop,rhot)
@ -1305,36 +1315,37 @@ contains
subroutine unset_eqspl
subroutine unset_equil_spline
implicit none
if(allocated(tr)) deallocate(tr)
if(allocated(tz)) deallocate(tz)
if(allocated(tfp)) deallocate(tfp)
if(allocated(cfp)) deallocate(cfp)
if(allocated(cceq)) deallocate(cceq)
if(allocated(cceq01)) deallocate(cceq01)
if(allocated(cceq10)) deallocate(cceq10)
if(allocated(cceq02)) deallocate(cceq02)
if(allocated(cceq20)) deallocate(cceq20)
if(allocated(cceq11)) deallocate(cceq11)
nsr=0
nsz=0
nsf=0
end subroutine unset_eqspl
if (allocated(tr)) deallocate(tr)
if (allocated(tz)) deallocate(tz)
if (allocated(tfp)) deallocate(tfp)
if (allocated(cfp)) deallocate(cfp)
if (allocated(cceq)) deallocate(cceq)
if (allocated(cceq01)) deallocate(cceq01)
if (allocated(cceq10)) deallocate(cceq10)
if (allocated(cceq02)) deallocate(cceq02)
if (allocated(cceq20)) deallocate(cceq20)
if (allocated(cceq11)) deallocate(cceq11)
nsr = 0
nsz = 0
nsf = 0
end subroutine unset_equil_spline
subroutine unset_q
implicit none
if(allocated(psinr)) deallocate(psinr)
if(allocated(cq)) deallocate(cq)
nq=0
if (allocated(psinr)) deallocate(psinr)
if (allocated(cq)) deallocate(cq)
nq = 0
end subroutine unset_q
subroutine unset_rhospl
subroutine unset_rho_spline
implicit none
if(allocated(rhopr)) deallocate(rhopr)
@ -1342,6 +1353,6 @@ contains
if(allocated(crhop)) deallocate(crhop)
if(allocated(crhot)) deallocate(crhot)
nrho=0
end subroutine unset_rhospl
end subroutine unset_rho_spline
end module equilibrium

18
src/gray_jetto1beam.f90
View File

@ -64,7 +64,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
! Set MHD equilibrium data
init_equilibrium: block
use equilibrium, only : set_eqspl
use equilibrium, only : set_equil_spline
! Copy argument arrays
data%equilibrium%rv = r
@ -81,12 +81,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
data%equilibrium%qpsi = qpsi
! Compute splines
call set_eqspl(params%equilibrium, data%equilibrium)
call set_equil_spline(params%equilibrium, data%equilibrium)
end block init_equilibrium
! Set plasma kinetic profiles
init_profiles: block
use coreprofiles, only : set_prfspl
use coreprofiles, only : set_profiles_spline
! Copy argument arrays
data%profiles%derad = dne
@ -94,7 +94,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
data%profiles%zfc = zeff
! Compute splines
call set_prfspl(params%profiles, data%profiles)
call set_profiles_spline(params%profiles, data%profiles)
end block init_profiles
! Set wave launcher parameters
@ -115,16 +115,16 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
! Free memory
free_memory: block
use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
use coreprofiles, only : unset_prfspl
use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q
use coreprofiles, only : unset_profiles_spline
! Unset global variables of the `equilibrium` module
call unset_eqspl
call unset_rhospl
call unset_equil_spline
call unset_rho_spline
call unset_q
! Unset global variables of the `coreprofiles` module
call unset_prfspl
call unset_profiles_spline
end block free_memory
! Copy over the results

57
src/main.f90
View File

@ -182,7 +182,7 @@ contains
! or an analytical description) and initialises the respective
! GRAY parameters and data.
use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
set_equian, set_eqspl, eq_scal
set_equil_an, set_equil_spline, eq_scal
implicit none
@ -190,17 +190,11 @@ contains
type(gray_parameters), intent(inout) :: params
type(gray_data), intent(out) :: data
if(params%equilibrium%iequil < 2) then
if (params%equilibrium%iequil < 2) then
! Analytical equilibrium
! TODO: rewrite using derived type
call read_equil_an(params%equilibrium%filenm, &
params%raytracing%ipass, &
data%equilibrium%rv, &
data%equilibrium%zv, &
data%equilibrium%fpol, &
data%equilibrium%qpsi, &
data%equilibrium%rlim, &
data%equilibrium%zlim)
data%equilibrium)
! Set psia sign to give the correct sign to Iphi
! (COCOS=3: psia<0 for Iphi>0)
@ -216,17 +210,10 @@ contains
call eq_scal(params%equilibrium, data%equilibrium)
! Set global variables (for splines)
if(params%equilibrium%iequil < 2) then
! TODO: rewrite using derived type
call set_equian(data%equilibrium%rv(1), &
data%equilibrium%zv(1), &
data%equilibrium%rv(2), &
data%equilibrium%fpol(1) / data%equilibrium%rv(1), &
data%equilibrium%qpsi(1), &
data%equilibrium%qpsi(2), &
data%equilibrium%qpsi(3))
if (params%equilibrium%iequil < 2) then
call set_equil_an(data%equilibrium)
else
call set_eqspl(params%equilibrium, data%equilibrium)
call set_equil_spline(params%equilibrium, data%equilibrium)
end if
end subroutine init_equilibrium
@ -234,7 +221,7 @@ contains
subroutine deinit_equilibrium(data)
! Free all memory allocated by the init_equilibrium subroutine.
use gray_params, only : equilibrium_data
use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q
implicit none
@ -248,8 +235,8 @@ contains
if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
! Unset global variables of the `equilibrium` module
call unset_eqspl
call unset_rhospl
call unset_equil_spline
call unset_rho_spline
call unset_q
end subroutine deinit_equilibrium
@ -261,7 +248,8 @@ contains
use gray_params, only : profiles_parameters, profiles_data
use equilibrium, only : frhopolv
use coreprofiles, only : read_profiles_an, read_profiles, &
scale_profiles, set_prfan, set_prfspl
scale_profiles, set_profiles_an, &
set_profiles_spline
implicit none
! subroutine arguments
@ -269,10 +257,9 @@ contains
real(wp_), intent(in) :: factb
type(profiles_data), intent(out) :: data
if(params%iprof == 0) then
if (params%iprof == 0) then
! Analytical profiles
! TODO: rewrite using derived type
call read_profiles_an(params%filenm, data%terad, data%derad, data%zfc)
call read_profiles_an(params%filenm, data)
else
! Numerical profiles
call read_profiles(params%filenm, data)
@ -291,13 +278,12 @@ contains
call scale_profiles(params, factb, data)
! Set global variables
if(params%iprof == 0) then
if (params%iprof == 0) then
! Analytical profiles
! TODO: rewrite using derived type
call set_prfan(data%terad, data%derad, data%zfc)
call set_profiles_an(data)
else
! Numerical profiles
call set_prfspl(params, data)
call set_profiles_spline(params, data)
end if
end subroutine init_profiles
@ -305,7 +291,7 @@ contains
subroutine deinit_profiles(data)
! Free all memory allocated by the init_profiles subroutine.
use gray_params, only : profiles_data
use coreprofiles, only : unset_prfspl
use coreprofiles, only : unset_profiles_spline
implicit none
@ -317,7 +303,7 @@ contains
if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
! Unset global variables of the `coreprofiles` module
call unset_prfspl
call unset_profiles_spline
end subroutine deinit_profiles
@ -375,7 +361,7 @@ contains
allocate(data%equilibrium%rlim(5))
allocate(data%equilibrium%zlim(5))
params%raytracing%ipass = abs(params%raytracing%ipass)
if(params%equilibrium%iequil < 2) then
if (params%equilibrium%iequil < 2) then
rmxm = (data%equilibrium%rv(1) + data%equilibrium%rv(2)) * 0.01_wp_
else
rmxm = data%equilibrium%rv(size(data%equilibrium%rv))
@ -406,7 +392,8 @@ contains
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx)
use const_and_precisions, only : zero, degree
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
use coreprofiles, only : set_profiles_an, set_profiles_spline, &
temp, fzeff
use dispersion, only : expinit
use gray_params, only : gray_parameters, print_parameters
use beams, only : launchangles2n, xgygcoeff
@ -457,7 +444,7 @@ contains
p0jk, ext, eyt, iiv)
! Initialise the dispersion module
if(params%ecrh_cd%iwarm > 1) call expinit
if (params%ecrh_cd%iwarm > 1) call expinit
! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot