diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index 024b9c4..6238b92 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -1,91 +1,149 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, & - p0mw, alphain, betain, dpdv, jcd, pabs, icd, ierr) - use const_and_precisions, only : wp_ - use units, only : ucenr,usumm,uprj0,uwbm,udisp,ubres,ucnt,uoutr,ueq,uprfin, & - uflx,upec,uprm,ubeam - use gray_params, only : read_params,raytracing,ecrh_cd,antenna,& - equilibrium,profiles,output - use beams, only : read_beam2 - use graycore, only : gray_main - use reflections, only : range2rect - use coreprofiles, only : tene_scal + p0mw, alphain, betain, dpdv, jcd, pabs, icd, error) + use const_and_precisions, only: wp_ + use gray_params, only: gray_parameters, gray_data, gray_results + use graycore, only: gray_main + implicit none -! arguments + + ! subroutine arguments integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam - real(wp_), dimension(mr), intent(in) :: r - real(wp_), dimension(mz), intent(in) :: z - real(wp_), dimension(mr,mz), intent(in) :: psin - real(wp_), intent(in) :: psia, rax, zax, p0mw, alphain, betain - real(wp_), dimension(nbnd), intent(in) :: rbnd, zbnd - real(wp_), dimension(nrho), intent(in) :: psrad, fpol, te, dne, zeff, qpsi - real(wp_), dimension(nrho), intent(out) :: dpdv, jcd - real(wp_), intent(out) :: pabs, icd - integer, intent(out) :: ierr -! local variables - real(wp_), dimension(nrho) :: psinr - integer :: iox0 - real(wp_) :: r0m,rvac,alpha0,beta0,psipol0,chipol0,rwallm - real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir - real(wp_), dimension(5) :: rlim,zlim - logical, save :: firstcall=.true. - type(raytracing) :: rtrp - type(ecrh_cd) :: hcdp - type(antenna) :: antp - type(equilibrium) :: eqp - type(profiles) :: prfp - type(output) :: outp + real(wp_), dimension(mr), intent(in) :: r + real(wp_), dimension(mz), intent(in) :: z + real(wp_), dimension(mr,mz), intent(in) :: psin + real(wp_), intent(in) :: psia, rax, zax, p0mw, alphain, betain + real(wp_), dimension(nbnd), intent(in) :: rbnd, zbnd + real(wp_), dimension(nrho), intent(in) :: psrad, fpol, te, dne, zeff, qpsi + real(wp_), dimension(nrho), intent(out) :: dpdv, jcd + real(wp_), intent(out) :: pabs, icd + integer, intent(out) :: error -! if first call read parameters from external file - if (firstcall) then - call read_params('gray.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp,uprm) - rwallm=r(1) - antp%filenm='graybeam.data' - eqp%filenm='JETTO' - eqp%iequil=ijetto+1 - prfp%filenm='JETTO' - outp%ipec=1 - firstcall=.false. - outp%nrho=nrho - end if - -! call tene_scal(te,dne,prfp%factte,prfp%factne,eqp%factb,& -! prfp%iscal,prfp%iprof) - - rvac=rax - psinr=psrad + ! local variables + type(gray_parameters) :: params + type(gray_data) :: data + type(gray_results) :: res + logical, save :: firstcall = .true. - alpha0=alphain - beta0=betain - call read_beam2(antp%filenm,ibeam,alpha0,beta0,fghz,antp%iox,x0,y0,z0, & - w1,w2,ri1,ri2,phiw,phir,ubeam) - psipol0=antp%psi - chipol0=antp%chi - iox0=antp%iox + ! Initialisation tasks for the first call only + first_call: if (firstcall) then + firstcall = .false. -! set simple limiter - r0m=sqrt(x0**2+y0**2)*0.01_wp_ - call range2rect(rwallm,max(r0m,r(mr)),z(1),z(mz),rlim,zlim) + ! Read parameters from external file + init_params: block + use gray_params, only : read_parameters, set_globals + call read_parameters('gray.data', params) -! call main subroutine for the ibeam-th beam - call gray_main(r,z,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & - psrad,te,dne,zeff,prfp, rlim,zlim, p0mw,fghz,alpha0,beta0,(/x0,y0,z0/), & - w1,w2,ri1,ri2,phiw,phir,iox0,psipol0,chipol0, dpdv,jcd,pabs,icd,outp, & - rtrp,hcdp,ierr,sqrt(psrad)) + ! Override some parameters + params%misc%rwall = r(1) + params%antenna%filenm = 'graybeam.data' + params%equilibrium%filenm = 'JETTO' + params%equilibrium%iequil = ijetto + 1 + params%profiles%filenm = 'JETTO' + params%output%ipec = 1 + params%output%nrho = nrho -! close output (debug) files - close(ucenr) - close(usumm) - close(uprj0) - close(uprj0+1) - close(uwbm) - close(udisp) - close(ubres) - close(ucnt) - close(uoutr) - close(ueq) - close(uprfin) - close(uflx) - close(upec) + ! Set the global variables of `gray_params` + call set_globals(params) + end block init_params + + ! Set a simple limiter + simple_limiter: block + use reflections, only : range2rect + use limiter, only : limiter_set_globals=>set_globals + + real(wp_) :: r0m + r0m = sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2)*0.01_wp_ + call range2rect(params%misc%rwall, max(r0m, r(mr)), z(1), z(mz), & + data%equilibrium%rlim, data%equilibrium%zlim) + ! Set the global variables of `limiter` + call limiter_set_globals(data%equilibrium) + end block simple_limiter + + end if first_call + + ! Set MHD equilibrium data + init_equilibrium: block + use equilibrium, only : set_eqspl + + ! Copy argument arrays + data%equilibrium%rv = r + data%equilibrium%zv = z + data%equilibrium%rax = rax + data%equilibrium%rvac = rax + data%equilibrium%zax = zax + data%equilibrium%psinr = psrad + data%equilibrium%fpol = fpol + data%equilibrium%rbnd = rbnd + data%equilibrium%zbnd = zbnd + data%equilibrium%psia = psia + data%equilibrium%psin = psin + data%equilibrium%qpsi = qpsi + + ! Compute splines + call set_eqspl(params%equilibrium, data%equilibrium) + end block init_equilibrium + + ! Set plasma kinetic profiles + init_profiles: block + use coreprofiles, only : set_prfspl + + ! Copy argument arrays + data%profiles%derad = dne + data%profiles%terad = te + data%profiles%zfc = zeff + + ! Compute splines + call set_prfspl(params%profiles, data%profiles) + end block init_profiles + + ! Set wave launcher parameters + init_antenna: block + use beams, only: read_beam2 + + ! Copy argument variables + params%antenna%alpha = alphain + params%antenna%beta = betain + params%antenna%power = p0mw + + ! Read beam description file + call read_beam2(params%antenna, beamid=ibeam) + end block init_antenna + + ! Call main subroutine for the ibeam-th beam + call gray_main(params, data, res, error, rhout=sqrt(psrad)) + + ! Free memory + free_memory: block + use equilibrium, only : unset_eqspl, unset_rhospl, unset_q + use coreprofiles, only : unset_prfspl + + ! Unset global variables of the `equilibrium` module + call unset_eqspl + call unset_rhospl + call unset_q + + ! Unset global variables of the `coreprofiles` module + call unset_prfspl + end block free_memory + + ! Close output units used for debugging + close_debug: block + use units, only : ucenr, usumm, uprj0, uwbm, udisp, ubres, & + ucnt, uoutr, ueq, uprfin, uflx, upec + integer :: i, debug_units(13) + + debug_units = [ucenr, usumm, uprj0, uprj0+1, uwbm, udisp, & + ubres, ucnt, uoutr, ueq, uprfin, uflx, upec] + do i=1,size(debug_units) + close(debug_units(i)) + end do + end block close_debug + + ! Copy over the results + pabs = res%pabs + icd = res%icd + dpdv = res%dpdv + jcd = res%jcd end subroutine gray_jetto1beam