From 948a512254cfbd4d576c4f4ec442db5089187681 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 15 Dec 2021 02:31:09 +0100 Subject: [PATCH] src: use derived type arguments (work in progress) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This change structures the arguments of most functions, in particular gray_main, into well-defined categories using derived types. All types are defined in the gray_params.f90 (location subject to change) and are organised as follows: gray_parameters (statically allocated data) ├── antenna_parameters ├── ecrh_cd_parameters ├── equilibrium_parameters ├── misc_parameters ├── output_parameters ├── profiles_parameters └── raytracing_parameters gray_data - inputs of gray_main (dynamically-allocated arrays) ├── equilibrium_data └── profiles_data gray_results - outputs of gray_main (dynamically-allocated arrays) --- src/beamdata.f90 | 11 +- src/beams.f90 | 487 +++++++++++++++--------------- src/equilibrium.f90 | 568 ++++++++++++++++++----------------- src/gray_jetto1beam.f90 | 16 +- src/gray_params.f90 | 395 ++++++++++++++----------- src/graycore.f90 | 542 +++++++++------------------------- src/limiter.f90 | 51 ++-- src/magsurf_data.f90 | 6 +- src/main.f90 | 638 +++++++++++++++++++++++++++------------- 9 files changed, 1400 insertions(+), 1314 deletions(-) diff --git a/src/beamdata.f90 b/src/beamdata.f90 index dc76472..298c26d 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -9,10 +9,10 @@ contains subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - use gray_params, only : rtrparam_type + use gray_params, only : raytracing_parameters use const_and_precisions, only : zero,half,two implicit none - type(rtrparam_type), intent(in) :: rtrparam + type(raytracing_parameters), intent(in) :: rtrparam real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & gri,psjki,ppabs,ccci real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri @@ -50,6 +50,10 @@ contains nstep=rtrparam%nstep + ! Allocate for each ray: + ! y = (x, k), + ! yp = dy/dt, + ! etc. call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) end subroutine init_btr @@ -97,6 +101,9 @@ contains + ! Functions to map radial (j), poloidal (k) ray + ! indices to a single global index (i) + function rayi2jk(i) result(jk) implicit none integer, intent(in) :: i diff --git a/src/beams.f90 b/src/beams.f90 index 2e0a247..25d988b 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -4,166 +4,175 @@ module beams contains - subroutine read_beam0(file_beam,fghz,x00,y00,z00, & - wcsi,weta,rcicsi,rcieta,phiw,phir,unit) - use const_and_precisions, only : pi,vc=>ccgs_ - use utils, only : get_free_unit + subroutine read_beam0(params, unit) + ! Reads the wave launcher parameters for the simple case + ! where w(z) and 1/R(z) are fixed. + + use const_and_precisions, only : pi, vc=>ccgs_ + use gray_params, only : antenna_parameters + use utils, only : get_free_unit + implicit none -! arguments - character(len=*), intent(in) :: file_beam - real(wp_), intent(out) :: fGHz,x00,y00,z00 - real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw + + ! subroutine arguments + type(antenna_parameters), intent(inout) :: params integer, intent(in), optional :: unit -! local variables + + ! local variables integer :: u - real(wp_) :: ak0,zrcsi,zreta,w0csi,w0eta,d0csi,d0eta - + real(wp_) :: ak0,zrcsi,zreta + if (present(unit)) then - u=unit + u = unit else u = get_free_unit() end if - open(unit=u,file=trim(file_beam),status='OLD',action='READ') -! fghz wave frequency (GHz) - read(u,*) fGHz -! x00,y00,z00 coordinates of launching point in cm - read(u,*) x00, y00, z00 -! beams parameters in local reference system -! w0 -> waist, d0 -> waist distance from launching point -! phiw angle of spot ellipse - read(u,*) w0csi,w0eta,d0csi,d0eta,phiw + open(unit=u, file=trim(params%filenm), status='OLD', action='READ') + read(u, *) params%fghz + read(u, *) params%pos + read(u, *) params%w, params%ri, params%phi(1) close(u) - ak0=2.0e9_wp_*pi*fghz/vc - zrcsi=0.5_wp_*ak0*w0csi**2 - zreta=0.5_wp_*ak0*w0eta**2 - - wcsi=w0csi*sqrt(1.0_wp_+(d0csi/zrcsi)**2) - weta=w0eta*sqrt(1.0_wp_+(d0eta/zreta)**2) - rcicsi=-d0csi/(d0csi**2+zrcsi**2) - rcieta=-d0eta/(d0eta**2+zreta**2) - phir=phiw + ak0 = 2.0e9_wp_* pi * params%fghz / vc + zrcsi = 0.5_wp_ * ak0 * params%w(1)**2 + zreta = 0.5_wp_ * ak0 * params%w(2)**2 + params%w(1) = params%w(1) * sqrt(1.0_wp_ + (params%ri(1)/zrcsi)**2) + params%w(2) = params%w(2) * sqrt(1.0_wp_ + (params%ri(2)/zreta)**2) + params%ri(1) = -params%ri(1) / (params%ri(1)**2 + zrcsi**2) + params%ri(2) = -params%ri(2) / (params%ri(2)**2 + zreta**2) + params%phi(2) = params%phi(1) end subroutine read_beam0 + subroutine read_beam1(params, unit) + ! Reads the wave launcher parameters for the case + ! where w(z, α) and 1/R(z, α) depend on the launcher angle α. - subroutine read_beam1(file_beam,alpha0,beta0,fghz,x00,y00,z00, & - wcsi,weta,rcicsi,rcieta,phiw,phir,unit) use const_and_precisions, only : pi,vc=>ccgs_ - use simplespline, only : spli, difcs - use utils, only : get_free_unit,locate + use gray_params, only : antenna_parameters + use simplespline, only : spli, difcs + use utils, only : get_free_unit,locate + implicit none -! arguments - character(len=*), intent(in) :: file_beam - real(wp_), intent(inout) :: alpha0 - real(wp_), intent(out) :: fghz,x00,y00,z00,beta0 - real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw + + ! subroutine arguments + type(antenna_parameters), intent(inout) :: params integer, intent(in), optional :: unit -! local variables - integer :: u,iopt,ier,nisteer,i,k,ii - real(wp_) :: steer,dal - real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & - z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & - cbeta,cx0,cy0,cz0,cwaist1,cwaist2, & - crci1,crci2,cphi1,cphi2 - + + ! local variables + integer :: u, iopt, ier, nisteer, i, k, ii + real(wp_) :: steer, dal + real(wp_), dimension(:), allocatable :: & + alphastv, betastv, x00v, y00v, & + z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, & + cbeta, cx0, cy0, cz0, cwaist1, cwaist2, & + crci1, crci2, cphi1, cphi2 + if (present(unit)) then - u=unit + u = unit else u = get_free_unit() end if - open(unit=u,file=file_beam,status='OLD',action='READ') - read(u,*) fghz - + open(unit=u, file=params%filenm, status='OLD', action='READ') + read(u,*) params%fghz read(u,*) nisteer - allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & - waist2v(nisteer),rci1v(nisteer),rci2v(nisteer), & - phi1v(nisteer),phi2v(nisteer),x00v(nisteer), & - y00v(nisteer),z00v(nisteer),cbeta(4*nisteer), & - cx0(4*nisteer),cy0(4*nisteer),cz0(4*nisteer), & - cwaist1(4*nisteer),cwaist2(4*nisteer),crci1(4*nisteer), & - crci2(4*nisteer),cphi1(4*nisteer),cphi2(4*nisteer)) + allocate(alphastv(nisteer), betastv(nisteer), waist1v(nisteer), & + waist2v(nisteer), rci1v(nisteer), rci2v(nisteer), & + phi1v(nisteer), phi2v(nisteer), x00v(nisteer), & + y00v(nisteer), z00v(nisteer), cbeta(4*nisteer), & + cx0(4*nisteer), cy0(4*nisteer), cz0(4*nisteer), & + cwaist1(4*nisteer), cwaist2(4*nisteer), crci1(4*nisteer), & + crci2(4*nisteer), cphi1(4*nisteer), cphi2(4*nisteer)) do i=1,nisteer - read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), & - waist1v(i),waist2v(i),rci1v(i),rci2v(i),phi1v(i),phi2v(i) + read(u, *) steer, alphastv(i), betastv(i), & + x00v(i), y00v(i), z00v(i), & + waist1v(i), waist2v(i), & + rci1v(i), rci2v(i), & + phi1v(i), phi2v(i) end do close(u) -! initial beam data measured in mm -> transformed to cm - x00v = 0.1_wp_*x00v - y00v = 0.1_wp_*y00v - z00v = 0.1_wp_*z00v - waist1v = 0.1_wp_*waist1v - waist2v = 0.1_wp_*waist2v - rci1v = 10._wp_*rci1v - rci2v = 10._wp_*rci2v - - iopt=0 - call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier) - call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier) - call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier) - call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier) - call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier) - call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier) - call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier) - call difcs(alphastv,x00v,nisteer,iopt,cx0,ier) - call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) - call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) - if((alpha0 > alphastv(1)).and.(alpha0 < alphastv(nisteer))) then - call locate(alphastv,nisteer,alpha0,k) - dal=alpha0-alphastv(k) - beta0=spli(cbeta,nisteer,k,dal) - x00=spli(cx0,nisteer,k,dal) - y00=spli(cy0,nisteer,k,dal) - z00=spli(cz0,nisteer,k,dal) - wcsi=spli(cwaist1,nisteer,k,dal) - weta=spli(cwaist2,nisteer,k,dal) - rcicsi=spli(crci1,nisteer,k,dal) - rcieta=spli(crci2,nisteer,k,dal) - phiw=spli(cphi1,nisteer,k,dal) - phir=spli(cphi2,nisteer,k,dal) + ! initial beam data measured in mm -> transformed to cm + x00v = 0.1_wp_ * x00v + y00v = 0.1_wp_ * y00v + z00v = 0.1_wp_ * z00v + waist1v = 0.1_wp_ * waist1v + waist2v = 0.1_wp_ * waist2v + rci1v = 10._wp_ * rci1v + rci2v = 10._wp_ * rci2v + + iopt = 0 + call difcs(alphastv, betastv, nisteer, iopt, cbeta, ier) + call difcs(alphastv, waist1v, nisteer, iopt, cwaist1, ier) + call difcs(alphastv, rci1v, nisteer, iopt, crci1, ier) + call difcs(alphastv, waist2v, nisteer, iopt, cwaist2, ier) + call difcs(alphastv, rci2v, nisteer, iopt, crci2, ier) + call difcs(alphastv, phi1v, nisteer, iopt, cphi1, ier) + call difcs(alphastv, phi2v, nisteer, iopt, cphi2, ier) + call difcs(alphastv, x00v, nisteer, iopt, cx0, ier) + call difcs(alphastv, y00v, nisteer, iopt, cy0, ier) + call difcs(alphastv, z00v, nisteer, iopt, cz0, ier) + + if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then + call locate(alphastv, nisteer, params%alpha , k) + dal = params%alpha - alphastv(k) + params%beta = spli(cbeta, nisteer, k, dal) + params%pos(1) = spli(cx0, nisteer, k, dal) + params%pos(2) = spli(cy0, nisteer, k, dal) + params%pos(3) = spli(cz0, nisteer, k, dal) + params%w(1) = spli(cwaist1, nisteer, k, dal) + params%w(2) = spli(cwaist2, nisteer, k, dal) + params%ri(1) = spli(crci1, nisteer, k, dal) + params%ri(2) = spli(crci2, nisteer, k, dal) + params%phi(1) = spli(cphi1, nisteer, k, dal) + params%phi(2) = spli(cphi2, nisteer, k, dal) else - ! alpha0 outside table range - if(alpha0 >= alphastv(nisteer)) ii=nisteer - if(alpha0 <= alphastv(1)) ii=1 - alpha0=alphastv(ii) - beta0=betastv(ii) - x00=x00v(ii) - y00=y00v(ii) - z00=z00v(ii) - wcsi=waist1v(ii) - weta=waist2v(ii) - rcicsi=rci1v(ii) - rcieta=rci2v(ii) - phiw=phi1v(ii) - phir=phi2v(ii) + ! params%alpha outside table range + if(params%alpha >= alphastv(nisteer)) ii=nisteer + if(params%alpha <= alphastv(1)) ii=1 + params%alpha = alphastv(ii) + params%beta = betastv(ii) + params%pos(1) = x00v(ii) + params%pos(2) = y00v(ii) + params%pos(3) = z00v(ii) + params%w(1) = waist1v(ii) + params%w(2) = waist2v(ii) + params%ri(1) = rci1v(ii) + params%ri(2) = rci2v(ii) + params%phi(1) = phi1v(ii) + params%phi(2) = phi2v(ii) end if - deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & - phi1v,phi2v,x00v,y00v,z00v,cbeta, & - cx0,cy0,cz0,cwaist1,cwaist2,crci1,crci2,cphi1,cphi2) + deallocate(alphastv, betastv, waist1v, waist2v, rci1v, rci2v, & + phi1v, phi2v, x00v, y00v, z00v, cbeta, & + cx0, cy0, cz0, cwaist1, cwaist2, & + crci1, crci2, cphi1, cphi2) end subroutine read_beam1 - subroutine read_beam2(file_beam,beamid,alpha0,beta0,fghz,iox,x00,y00,z00, & - wcsi,weta,rcicsi,rcieta,phiw,phir,unit) - use utils, only : get_free_unit, intlin, locate + subroutine read_beam2(params, beamid, unit) + ! Reads the wave launcher parameters for the general case + ! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β. + + use gray_params, only : antenna_parameters + use utils, only : get_free_unit, intlin, locate use reflections, only : inside - use dierckx, only : curfit, splev, surfit, bispev + use dierckx, only : curfit, splev, surfit, bispev + implicit none - character(len=*), intent(in) :: file_beam + + ! subroutine arguments + type(antenna_parameters), intent(inout) :: params integer, intent(in) :: beamid - real(wp_), intent(inout) :: alpha0,beta0 - real(wp_), intent(out) :: fghz,x00,y00,z00, wcsi,weta,rcicsi,rcieta,phir,phiw - integer, intent(out) :: iox integer, intent(in), optional :: unit + ! local variables character(len=20) :: beamname integer :: u integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta @@ -189,11 +198,12 @@ contains real(wp_), parameter :: sspl=0.01_wp_ if (present(unit)) then - u=unit + u = unit else u = get_free_unit() end if - open(unit=u,file=file_beam,status='OLD',action='READ') + + open(unit=u, file=params%filenm, status='OLD', action='READ') !======================================================================================= ! # of beams read(u,*) nbeam @@ -202,13 +212,13 @@ contains jumprow=0 ! c==================================================================================== do i=1,beamid-1 - read(u,*) beamname, iox, fghz, nalpha, nbeta + read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta jumprow = jumprow+nalpha*nbeta end do ! c==================================================================================== ! ! beam of interest - read(u,*) beamname, iox, fghz, nalpha, nbeta + read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta ! ! c==================================================================================== ! unused beams' data grids @@ -216,7 +226,7 @@ contains read(u,*) beamname end do do i=1,jumprow - read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 + read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2 end do ! c==================================================================================== ! @@ -231,12 +241,12 @@ contains ! c==================================================================================== ! beam data grid reading do i=1,nisteer - read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 + read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2 ! -! initial beam data (x00, y00, z00) are measured in mm -> transformed to cm - x00v(i)=0.1d0*x00 - y00v(i)=0.1d0*y00 - z00v(i)=0.1d0*z00 +! initial beam data (params%pos(1), params%pos(2), params%pos(3)) are measured in mm -> transformed to cm + x00v(i)=0.1d0*params%pos(1) + y00v(i)=0.1d0*params%pos(2) + z00v(i)=0.1d0*params%pos(3) alphastv(i)=alphast betastv(i)=betast waist1v(i)=0.1d0*waist1 @@ -259,17 +269,17 @@ contains ! ! no free variables if(fdeg.eq.3) then - alpha0=alphastv(1) - beta0=betastv(1) - x00=x00v(1) - y00=y00v(1) - z00=z00v(1) - wcsi=waist1v(1) - weta=waist2v(1) - rcicsi=rci1v(1) - rcieta=rci2v(1) - phiw=phi1v(1) - phir=phi2v(1) + params%alpha=alphastv(1) + params%beta=betastv(1) + params%pos(1)=x00v(1) + params%pos(2)=y00v(1) + params%pos(3)=z00v(1) + params%w(1)=waist1v(1) + params%w(2)=waist2v(1) + params%ri(1)=rci1v(1) + params%ri(2)=rci2v(1) + params%phi(2)=phi1v(1) + params%phi(1)=phi2v(1) deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & phi2v,x00v,y00v,z00v,xcoord,ycoord) return @@ -283,8 +293,8 @@ contains ! alpha = dependent variable xcoord = betastv ycoord = alphastv - xcoord0 = beta0 - ycoord0 = alpha0 + xcoord0 = params%beta + ycoord0 = params%alpha kx=min(nbeta-1,kspl) ! c==================================================================================== else @@ -293,8 +303,8 @@ contains ! beta = dependent/independent (fdeg = 1/0) xcoord = alphastv ycoord = betastv - xcoord0 = alpha0 - ycoord0 = beta0 + xcoord0 = params%alpha + ycoord0 = params%beta nxcoord = nalpha nycoord = nbeta kx=min(nalpha-1,kspl) @@ -461,23 +471,23 @@ contains call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier) ycoord0=fi(1) call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier) - wcsi=fi(1) + params%w(1)=fi(1) call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier) - weta=fi(1) + params%w(2)=fi(1) call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier) - rcicsi=fi(1) + params%ri(1)=fi(1) call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier) - rcieta=fi(1) + params%ri(2)=fi(1) call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier) - phiw=fi(1) + params%phi(2)=fi(1) call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier) - phir=fi(1) + params%phi(1)=fi(1) call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier) - x00=fi(1) + params%pos(1)=fi(1) call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier) - y00=fi(1) + params%pos(2)=fi(1) call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier) - z00=fi(1) + params%pos(3)=fi(1) ! c---------------------------------------------------------------------------------- else ! c---------------------------------------------------------------------------------- @@ -486,15 +496,15 @@ contains ! xcoord0=xcoord(ii) ycoord0=ycoord(ii) - x00=x00v(ii) - y00=y00v(ii) - z00=z00v(ii) - wcsi=waist1v(ii) - weta=waist2v(ii) - rcicsi=rci1v(ii) - rcieta=rci2v(ii) - phiw=phi1v(ii) - phir=phi2v(ii) + params%pos(1)=x00v(ii) + params%pos(2)=y00v(ii) + params%pos(3)=z00v(ii) + params%w(1)=waist1v(ii) + params%w(2)=waist2v(ii) + params%ri(1)=rci1v(ii) + params%ri(2)=rci2v(ii) + params%phi(2)=phi1v(ii) + params%phi(1)=phi2v(ii) end if ! c==================================================================================== else @@ -586,53 +596,53 @@ contains ! 5: xcoord0 unchanged, ycoord0 moved on side C ! 7: xcoord0 moved on side D, ycoord0 unchanged ! 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates -! in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in +! in 1,3,5,7 incheck is set back to 1 to evaluate params%pos(1),params%pos(2),params%pos(3),waist,rci,phi in ! new (xcoord0,ycoord0) -! in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the +! in 2,4,6,8 incheck remains 0 and params%pos(1),params%pos(2),params%pos(3),waist,rci,phi values at the ! (xcoord0,ycoord0) vertex are used - alpha0 = xcoord0 - beta0 = ycoord0 + params%alpha = xcoord0 + params%beta = ycoord0 SELECT CASE (in) CASE (1) - ! beta0 outside table range + ! params%beta outside table range ! locate position of xcoord0 with respect to x coordinates of side A call locate(xpolygA,nxcoord,xcoord0,ii) ! find corresponding y value on side A for xcoord position call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0) incheck = 1 CASE (2) - ! alpha0 and beta0 outside table range + ! params%alpha and params%beta outside table range ! xcoord0, ycoord0 set xcoord0 = xvert(2) ycoord0 = yvert(2) ii = nxcoord !indice per assegnare valori waist, rci, phi CASE (3) - ! alpha0 outside table range + ! params%alpha outside table range call locate(ypolygB,nycoord,ycoord0,ii) call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0) incheck = 1 CASE (4) - ! alpha0 and beta0 outside table range + ! params%alpha and params%beta outside table range xcoord0 = xvert(3) ycoord0 = yvert(3) ii = nxcoord+nycoord-1 CASE (5) - ! beta0 outside table range + ! params%beta outside table range call locate(xpolygC,nxcoord,xcoord0,ii) call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0) incheck = 1 CASE (6) - ! alpha0 and beta0 outside table range + ! params%alpha and params%beta outside table range xcoord0 = xvert(4) ycoord0 = yvert(4) ii = 2*nxcoord+nycoord-2 CASE (7) - ! alpha0 outside table range + ! params%alpha outside table range call locate(ypolygD,nycoord,ycoord0,ii) call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0) incheck = 1 CASE (8) - ! alpha0 and beta0 outside table range + ! params%alpha and params%beta outside table range xcoord0 = xvert(1) ycoord0 = yvert(1) ii = 1 @@ -651,44 +661,44 @@ contains allocate(wrk(lwrk),iwrk(kwrk)) call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - wcsi=fi(1) + params%w(1)=fi(1) call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - weta=fi(1) + params%w(2)=fi(1) call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - rcicsi=fi(1) + params%ri(1)=fi(1) call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - rcieta=fi(1) + params%ri(2)=fi(1) call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - phiw=fi(1) + params%phi(2)=fi(1) call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - phir=fi(1) + params%phi(1)=fi(1) call bispev(txx0,nxx0,tyx0,nyx0,cx0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - x00=fi(1) + params%pos(1)=fi(1) call bispev(txy0,nxy0,tyy0,nyy0,cy0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - y00=fi(1) + params%pos(2)=fi(1) call bispev(txz0,nxz0,tyz0,nyz0,cz0, & kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) - z00=fi(1) + params%pos(3)=fi(1) deallocate(wrk,iwrk) ! c---------------------------------------------------------------------------------- else ! c---------------------------------------------------------------------------------- - x00=x00v(ii) - y00=y00v(ii) - z00=z00v(ii) - wcsi=waist1v(ii) - weta=waist2v(ii) - rcicsi=rci1v(ii) - rcieta=rci2v(ii) - phiw=phi1v(ii) - phir=phi2v(ii) + params%pos(1)=x00v(ii) + params%pos(2)=y00v(ii) + params%pos(3)=z00v(ii) + params%w(1)=waist1v(ii) + params%w(2)=waist2v(ii) + params%ri(1)=rci1v(ii) + params%ri(2)=rci2v(ii) + params%phi(2)=phi1v(ii) + params%phi(1)=phi2v(ii) end if ! c==================================================================================== end if @@ -709,11 +719,11 @@ contains !####################################################################################### ! set correct values for alpha, beta if(fdeg.eq.2) then - alpha0 = ycoord0 - beta0 = xcoord0 + params%alpha = ycoord0 + params%beta = xcoord0 else - alpha0 = xcoord0 - beta0 = ycoord0 + params%alpha = xcoord0 + params%beta = ycoord0 end if !####################################################################################### deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & @@ -722,53 +732,62 @@ contains end subroutine read_beam2 - subroutine launchangles2n(alpha,beta,xv,anv) - use const_and_precisions, only : degree - implicit none -! arguments - real(wp_), intent(in) :: alpha,beta,xv(3) - real(wp_), intent(out) :: anv(3) -! local variables - real(wp_) :: r,anr,anphi,a,b + subroutine launchangles2n(params, anv) + ! Given the wave launcher `params` computes the initial + ! wavevector `anv`, defined by n̅ = ck̅/ω, in cartesian coordinates. - r=sqrt(xv(1)**2+xv(2)**2) -! phi=atan2(y,x) - a = degree*alpha - b = degree*beta -! -! angles alpha, beta in a local reference system as proposed by Gribov et al -! + use const_and_precisions, only : degree + use gray_params, only : antenna_parameters + + implicit none + + ! subroutine arguments + type(antenna_parameters), intent(in) :: params + real(wp_), intent(out) :: anv(3) + + ! local variables + real(wp_) :: r, anr, anphi, a, b + + r = sqrt(params%pos(1)**2 + params%pos(2)**2) + a = degree*params%alpha + b = degree*params%beta + + ! Angles α, β in a local reference system + ! as proposed by Gribov et al. anr = -cos(b)*cos(a) anphi = sin(b) -! anx = -cos(b)*cos(a) -! any = sin(b) - anv(1) = (anr*xv(1) - anphi*xv(2))/r ! = anx - anv(2) = (anr*xv(2) + anphi*xv(1))/r ! = any -! anr = (anx*xv(1) + any*xv(2))/r -! anphi = (any*xv(1) - anx*xv(2))/r - - anv(3) =-cos(b)*sin(a) ! = anz + anv(1) = (anr*params%pos(1) - anphi*params%pos(2))/r ! = anx + anv(2) = (anr*params%pos(2) + anphi*params%pos(1))/r ! = any + anv(3) = -cos(b)*sin(a) ! = anz end subroutine launchangles2n - subroutine xgygcoeff(fghz,ak0,bres,xgcn) - use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_,pi,wce1_ + + subroutine xgygcoeff(fghz, ak0, bres, xgcn) + ! Given the EC wave frequency computes: + ! + ! 1. vacuum wavevector `k0` (k₀ = ω/c), + ! 2. resonant magnetic field `bres` (qB/m = ω), + ! 3. adimensional `xgcn` parameter (X = ω_p²/ω² = nq²/ε₀mω²). + use const_and_precisions, only : qe=>ecgs_, me=>mecgs_, & + vc=>ccgs_, pi, wce1_ implicit none -! arguments - real(wp_), intent(in) :: fghz - real(wp_), intent(out) :: ak0,bres,xgcn -! local variables + + ! subroutine arguments + real(wp_), intent(in) :: fghz + real(wp_), intent(out) :: ak0, bres, xgcn + + ! local variables real(wp_) :: omega - omega=2.0e9_wp_*pi*fghz ! [rad/s] - ak0=omega/vc ! [rad/cm] -! -! yg=btot/bres -! - bres=omega/wce1_ ! [T] -! -! xg=xgcn*dens19 -! - xgcn=4.0e13_wp_*pi*qe**2/(me*omega**2) ! [10^-19 m^3] + omega = 2.0e9_wp_*pi*fghz ! [rad/s] + ak0 = omega/vc ! [rad/cm] + + ! yg = btot/bres + bres = omega/wce1_ ! [T] + + ! xg = xgcn*dens19 + xgcn = 4.0e13_wp_ * pi * qe**2/(me * omega**2) ! [10^-19 m^3] end subroutine xgygcoeff + end module beams diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 9a92494..1ee7085 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -13,7 +13,7 @@ module equilibrium ! === 2d spline psi(r,z), normalization and derivatives ========== integer, save :: nsr, nsz - real(wp_), save :: psia, psiant, psinop, psnbd + real(wp_), save :: psia, psiant, psinop, psnbnd real(wp_), dimension(:), allocatable, save :: tr,tz real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, & cceq20, cceq02, cceq11 @@ -34,121 +34,127 @@ module equilibrium contains - subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, & - rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit) + subroutine read_eqdsk(params, data, unit) + ! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm). + ! If given, the file is opened in the `unit` number. + ! For a description of the G-EQDSK, see the GRAY user manual. use const_and_precisions, only : one - use utils, only : get_free_unit + use gray_params, only : equilibrium_parameters, equilibrium_data + use utils, only : get_free_unit + implicit none -! arguments - character(len=*), intent(in) :: filenm - real(wp_), intent(out) :: psia,rvac,rax,zax - real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q - real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim - real(wp_), dimension(:,:), allocatable, intent(out) :: psin - integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit -! local variables - integer, parameter :: indef=0,iddef=1,iffdef=0 - integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim + + ! subroutine arguments + type(equilibrium_parameters), intent(in) :: params + type(equilibrium_data), intent(out) :: data + integer, optional, intent(in) :: unit + + ! local variables + integer :: u, idum, i, j, nr, nz, nbnd, nlim character(len=48) :: string - real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis - -! set default values if optional arguments are absent - in=indef; id=iddef; iffmt=iffdef - if(present(ipsinorm)) in=ipsinorm - if(present(idesc)) id=idesc - if(present(ifreefmt)) iffmt=ifreefmt - if (present(unit)) then - u=unit + real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis + real(wp_) :: xdum ! dummy variable, used to discard data + + if(present(unit)) then + u = unit else - u=get_free_unit() + u = get_free_unit() end if -! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html) - open(file=trim(filenm),status='old',action='read',unit=u) + ! Open the G-EQDSK file + open(file=trim(params%filenm), status='old', action='read', unit=u) -! get size of main arrays and allocate them - if (id==1) then + ! get size of main arrays and allocate them + if (params%idesc == 1) then read (u,'(a48,3i4)') string,idum,nr,nz else - read (u,*) nr,nz + read (u,*) nr, nz end if - if (allocated(rv)) deallocate(rv) - if (allocated(zv)) deallocate(zv) - if (allocated(psin)) deallocate(psin) - if (allocated(psinr)) deallocate(psinr) - if (allocated(fpol)) deallocate(fpol) - if (allocated(q)) deallocate(q) - allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr)) + if (allocated(data%rv)) deallocate(data%rv) + if (allocated(data%zv)) deallocate(data%zv) + if (allocated(data%psin)) deallocate(data%psin) + if (allocated(data%psinr)) deallocate(data%psinr) + if (allocated(data%fpol)) deallocate(data%fpol) + if (allocated(data%qpsi)) deallocate(data%qpsi) + allocate(data%rv(nr), data%zv(nz), & + data%psin(nr, nz), & + data%psinr(nr), & + data%fpol(nr), & + data%qpsi(nr)) -! store 0D data and main arrays - if (iffmt==1) then - read (u,*) dr,dz,rvac,rleft,zmid - read (u,*) rax,zax,psiaxis,psiedge,xdum - read (u,*) xdum,xdum,xdum,xdum,xdum - read (u,*) xdum,xdum,xdum,xdum,xdum - read (u,*) (fpol(i),i=1,nr) - read (u,*) (xdum,i=1,nr) - read (u,*) (xdum,i=1,nr) - read (u,*) (xdum,i=1,nr) - read (u,*) ((psin(i,j),i=1,nr),j=1,nz) - read (u,*) (q(i),i=1,nr) + ! Store 0D data and main arrays + if (params%ifreefmt==1) then + read (u, *) dr, dz, data%rvac, rleft, zmid + read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum + read (u, *) xdum, xdum, xdum, xdum, xdum + read (u, *) xdum, xdum, xdum, xdum, xdum + read (u, *) (data%fpol(i), i=1,nr) + read (u, *) (xdum,i=1, nr) + read (u, *) (xdum,i=1, nr) + read (u, *) (xdum,i=1, nr) + read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz) + read (u, *) (data%qpsi(i), i=1,nr) else - read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid - read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum - read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum - read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum - read (u,'(5e16.9)') (fpol(i),i=1,nr) - read (u,'(5e16.9)') (xdum,i=1,nr) - read (u,'(5e16.9)') (xdum,i=1,nr) - read (u,'(5e16.9)') (xdum,i=1,nr) - read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz) - read (u,'(5e16.9)') (q(i),i=1,nr) + read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid + read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum + read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum + read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum + read (u, '(5e16.9)') (data%fpol(i),i=1,nr) + read (u, '(5e16.9)') (xdum,i=1,nr) + read (u, '(5e16.9)') (xdum,i=1,nr) + read (u, '(5e16.9)') (xdum,i=1,nr) + read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz) + read (u, '(5e16.9)') (data%qpsi(i),i=1,nr) end if -! get size of boundary and limiter arrays and allocate them - read (u,*) nbnd,nlim - if (allocated(rbnd)) deallocate(rbnd) - if (allocated(zbnd)) deallocate(zbnd) - if (allocated(rlim)) deallocate(rlim) - if (allocated(zlim)) deallocate(zlim) + ! Get size of boundary and limiter arrays and allocate them + read (u,*) nbnd, nlim + if (allocated(data%rbnd)) deallocate(data%rbnd) + if (allocated(data%zbnd)) deallocate(data%zbnd) + if (allocated(data%rlim)) deallocate(data%rlim) + if (allocated(data%zlim)) deallocate(data%zlim) -! store boundary and limiter data - if(nbnd>0) then - allocate(rbnd(nbnd),zbnd(nbnd)) - if (iffmt==1) then - read(u,*) (rbnd(i),zbnd(i),i=1,nbnd) + ! Load plasma boundary data + if(nbnd > 0) then + allocate(data%rbnd(nbnd), data%zbnd(nbnd)) + if (params%ifreefmt == 1) then + read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd) else - read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd) - end if - end if - if(nlim>0) then - allocate(rlim(nlim),zlim(nlim)) - if (iffmt==1) then - read(u,*) (rlim(i),zlim(i),i=1,nlim) - else - read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim) + read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd) end if end if -! reading of G EQDSK file completed + ! Load limiter data + if(nlim > 0) then + allocate(data%rlim(nlim), data%zlim(nlim)) + if (params%ifreefmt == 1) then + read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim) + else + read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim) + end if + end if + + ! End of G-EQDSK file close(u) -! build rv,zv,psinr arrays and normalize psin - zleft=zmid-0.5_wp_*dz - dr=dr/(nr-1) - dz=dz/(nz-1) - dps=one/(nr-1) + ! Build rv,zv,psinr arrays + zleft = zmid-0.5_wp_*dz + dr = dr/(nr-1) + dz = dz/(nz-1) + dps = one/(nr-1) do i=1,nr - psinr(i)=(i-1)*dps - rv(i)=rleft+(i-1)*dr + data%psinr(i) = (i-1)*dps + data%rv(i) = rleft + (i-1)*dr end do do i=1,nz - zv(i)=zleft+(i-1)*dz + data%zv(i) = zleft + (i-1)*dz end do - psia=psiedge-psiaxis - if(in==0) psin=(psin-psiaxis)/psia - end subroutine read_eqdsk + ! Normalize psin + data%psia = psiedge - psiaxis + if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia + + end subroutine read_eqdsk subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit) @@ -211,191 +217,228 @@ contains close(u) end subroutine read_equil_an - subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr) - use const_and_precisions, only : zero,one,pi + + subroutine change_cocos(data, cocosin, cocosout, error) + ! Convert the MHD equilibrium data from one coordinate convention + ! (COCOS) to another. These are specified by `cocosin` and + ! `cocosout`, respectively. + ! + ! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010 + use const_and_precisions, only : zero, one, pi + use gray_params, only : equilibrium_data + implicit none -! arguments - real(wp_), intent(inout) :: psia - real(wp_), dimension(:), intent(inout) :: fpol,q + + ! subroutine arguments + type(equilibrium_data), intent(inout) :: data integer, intent(in) :: cocosin, cocosout - integer, intent(out), optional :: ierr -! local variables - real(wp_) :: isign,bsign - integer :: exp2pi,exp2piout - logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout + integer, intent(out), optional :: error - call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos) - call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout) + ! local variables + real(wp_) :: isign, bsign + integer :: exp2pi, exp2piout + logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout -! check sign consistency - isign=sign(one,psia) - if (.not.psiincr) isign=-isign - bsign=sign(one,fpol(size(fpol))) - if (qpos.neqv.isign*bsign*q(size(q))>zero) then -! warning: sign inconsistency found among q, Ipla and Bref - q=-q - if(present(ierr)) ierr=1 + call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos) + call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout) + + ! Check sign consistency + isign = sign(one, data%psia) + if (.not.psiincr) isign = -isign + bsign = sign(one, data%fpol(size(data%fpol))) + if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then + ! Warning: sign inconsistency found among q, Ipla and Bref + data%qpsi = -data%qpsi + if (present(error)) error = 1 else - if(present(ierr)) ierr=0 + if (present(error)) error = 0 end if -! convert cocosin to cocosout + ! Convert cocosin to cocosout -! opposite direction of toroidal angle phi in cocosin and cocosout - if (phiccw.neqv.phiccwout) fpol=-fpol -! q has opposite sign for given sign of Bphi*Ip - if (qpos .neqv. qposout) q=-q -! psi and Ip signs don't change accordingly - if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia -! convert Wb to Wb/rad or viceversa - psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi) + ! Opposite direction of toroidal angle phi in cocosin and cocosout + if (phiccw .neqv. phiccwout) data%fpol = -data%fpol + + ! q has opposite sign for given sign of Bphi*Ip + if (qpos .neqv. qposout) data%qpsi = -data%qpsi + + ! psi and Ip signs don't change accordingly + if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) & + data%psia = -data%psia + + ! Convert Wb to Wb/rad or viceversa + data%psia = data%psia * (2.0_wp_*pi)**(exp2piout - exp2pi) end subroutine change_cocos - subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos) - implicit none - integer, intent(in) :: cocos - integer, intent(out) :: exp2pi - logical, intent(out) :: phiccw,psiincr,qpos - integer :: cmod10,cmod4 - cmod10=mod(cocos,10) - cmod4=mod(cmod10,4) -! cocos>10 psi in Wb, cocos<10 psi in Wb/rad - exp2pi=(cocos-cmod10)/10 -! cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW - phiccw=(mod(cmod10,2)==1) -! cocos mod 10 = 1,2,5,6: psi increasing with positive Ip - psiincr=(cmod4==1 .or. cmod4==2) -! cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip - qpos=(cmod10<3 .or. cmod10>6) + subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos) + implicit none + + ! subroutine arguments + integer, intent(in) :: cocos + integer, intent(out) :: exp2pi + logical, intent(out) :: phiccw, psiincr, qpos + + ! local variables + integer :: cmod10, cmod4 + + cmod10 = mod(cocos, 10) + cmod4 = mod(cmod10, 4) + + ! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad + exp2pi = (cocos - cmod10)/10 + + ! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW + phiccw = (mod(cmod10, 2)== 1) + + ! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip + psiincr = (cmod4==1 .or. cmod4==2) + + ! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip + qpos = (cmod10<3 .or. cmod10>6) end subroutine decode_cocos - subroutine eq_scal(psia,fpol,isign,bsign,factb) + + subroutine eq_scal(params, data) + ! Rescale the magnetic field (B) and the plasma current (I) + ! and/or force their signs. + ! + ! Notes: + ! 1. signi and signb are ignored on input if equal to 0. + ! They are used to assign the direction of Bphi and Ipla BEFORE scaling. + ! 2. cocos=3 assumed: CCW direction is >0 + ! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged + ! 4. factb<0 reverses the directions of Bphi and Ipla + use const_and_precisions, only : one + use gray_params, only : equilibrium_parameters, equilibrium_data + implicit none - real(wp_), intent(inout) :: psia - integer, intent(inout) :: isign,bsign - real(wp_), dimension(:), intent(inout) :: fpol - real(wp_), intent(in) :: factb + + ! subroutine arguments + type(equilibrium_parameters), intent(inout) :: params + type(equilibrium_data), intent(inout) :: data + + ! local variables + real(wp_) :: last_fpol + + last_fpol = data%fpol(size(data%fpol)) - ! isign and bsign ignored on input if equal to 0 - ! used to assign the direction of Bphi and Ipla BEFORE scaling - ! cocos=3 assumed: CCW direction is >0 - ! Bphi and Ipla scaled by the same factor factb to keep q unchanged - ! factb<0 reverses the directions of Bphi and Ipla - if(isign/=0) psia=sign(psia,real(-isign,wp_)) - if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol - psia=psia*factb - fpol=fpol*factb - isign=int(sign(one,-psia)) - bsign=int(sign(one,fpol(size(fpol)))) + if (params%sgni /=0) & + data%psia = sign(data%psia, real(-params%sgni, wp_)) + if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) & + data%fpol = -data%fpol + + data%psia = data%psia * params%factb + data%fpol = data%fpol * params%factb + params%sgni = int(sign(one, -data%psia)) + params%sgnb = int(sign(one, last_fpol)) end subroutine eq_scal - subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,qpsi,sspl,ssfp, & - r0,rax,zax,rbnd,zbnd,ixp) - use const_and_precisions, only : zero,one - use dierckx, only : regrid,coeff_parder,curfit,splev - use gray_params, only : iequil - use reflections, only : inside - use utils, only : vmaxmin,vmaxmini + + subroutine set_eqspl(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 + use gray_params, only : equilibrium_parameters, equilibrium_data + use dierckx, only : regrid, coeff_parder, curfit, splev + use gray_params, only : iequil + use reflections, only : inside + use utils, only : vmaxmin, vmaxmini + implicit none -! arguments - real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol,qpsi - real(wp_), dimension(:,:), intent(in) :: psin - real(wp_), intent(in) :: psiwbrad - real(wp_), intent(in) :: sspl,ssfp - real(wp_), intent(in), optional :: r0,rax,zax - real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd - integer, intent(in), optional :: ixp -! local constants + + ! subroutine arguments + type(equilibrium_parameters), intent(in) :: params + type(equilibrium_data), intent(in) :: data + + ! local constants integer, parameter :: iopt=0 -! local variables + + ! local variables integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 - real(wp_), dimension(size(psinr)) :: rhotn + real(wp_), dimension(size(data%psinr)) :: rhotn real(wp_), dimension(1) :: fpoli real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk integer :: ier,ixploc,info,i,j,ij -! compute array sizes and prepare working space arrays - nr=size(rv) - nz=size(zv) + ! compute array sizes and prepare working space arrays + nr=size(data%rv) + nz=size(data%zv) nrz=nr*nz nrest=nr+ksplp nzest=nz+ksplp lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest) liwrk=nz+nr+nrest+nzest+kspl - npsi=size(psinr) + npsi=size(data%psinr) npsest=npsi+ksplp lwrkf=npsi*ksplp+npsest*(7+3*kspl) allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest))) -! spline fitting/interpolation of psin(i,j) and derivatives + ! spline fitting/interpolation of data%psin(i,j) and derivatives -! allocate knots and spline coefficients arrays + ! allocate knots and spline coefficients arrays if (allocated(tr)) deallocate(tr) if (allocated(tz)) deallocate(tz) if (allocated(cceq)) deallocate(cceq) allocate(tr(nrest),tz(nzest),cceq(nrz)) -! length in m !!! - - rmnm=rv(1) - rmxm=rv(nr) - zmnm=zv(1) - zmxm=zv(nz) + ! length in m !!! + rmnm=data%rv(1) + rmxm=data%rv(nr) + zmnm=data%zv(1) + zmxm=data%zv(nz) if (iequil>2) then -! data valid only inside boundary (psin=0 outside), e.g. source==ESCO -! presence of boundary anticipated here to filter invalid data - if(present(rbnd).and.present(zbnd)) then - nbnd=min(size(rbnd),size(zbnd)) - else - nbnd=0 - end if -! determine number of valid grid points + ! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO + ! presence of boundary anticipated here to filter invalid data + nbnd = min(size(data%rbnd), size(data%zbnd)) + + ! determine number of valid grid points nrz=0 do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle + if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle else - if(psin(i,j).le.0.0d0) cycle + if(data%psin(i,j).le.0.0d0) cycle end if nrz=nrz+1 end do end do -! store valid data + + ! store valid data allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) ij=0 do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle + if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle else - if(psin(i,j).le.0.0d0) cycle + if(data%psin(i,j).le.0.0d0) cycle end if ij=ij+1 - rv1d(ij)=rv(i) - zv1d(ij)=zv(j) - fvpsi(ij)=psin(i,j) + rv1d(ij)=data%rv(i) + zv1d(ij)=data%zv(j) + fvpsi(ij)=data%psin(i,j) wf(ij)=1.0d0 end do end do -! fit as a scattered set of points -! use reduced number of knots to limit memory comsumption ? + ! fit as a scattered set of points + ! use reduced number of knots to limit memory comsumption ? nsr=nr/4+4 nsz=nz/4+4 - sspln=sspl + sspln=params%ssplps call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) -! if ier=-1 data are fitted using sspl=0 + ! if ier=-1 data are fitted using params%ssplps=0 if(ier.eq.-1) then sspln=0.0_wp_ nsr=nr/4+4 @@ -404,29 +447,29 @@ contains rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) end if deallocate(rv1d,zv1d,wf,fvpsi) -! reset nrz to the total number of grid points for next allocations + ! reset nrz to the total number of grid points for next allocations nrz=nr*nz else -! iequil==2: data are valid on the full R,z grid + ! iequil==2: data are valid on the full R,z grid -! reshape 2D psi array to 1D (transposed) array and compute spline coeffs + ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs allocate(fvpsi(nrz)) - fvpsi=reshape(transpose(psin),(/nrz/)) - sspln=sspl - call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & + fvpsi=reshape(transpose(data%psin),(/nrz/)) + sspln=params%ssplps + call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & + kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) -! if ier=-1 data are re-fitted using sspl=0 + ! if ier=-1 data are re-fitted using params%ssplps=0 if(ier==-1) then sspln=0.0_wp_ - call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & + call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) end if deallocate(fvpsi) end if -! compute spline coefficients for psi partial derivatives + ! compute spline coefficients for psi partial derivatives lw10 = nr*(ksplp-1) + nz*ksplp + nrz lw01 = nr*ksplp + nz*(ksplp-1) + nrz lw20 = nr*(ksplp-2) + nz*ksplp + nrz @@ -444,74 +487,57 @@ contains call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier) -! spline interpolation of fpol(i) + ! Spline interpolation of data%fpol(i) -! allocate knots and spline coefficients arrays + ! Allocate knots and spline coefficients arrays if (allocated(tfp)) deallocate(tfp) if (allocated(cfp)) deallocate(cfp) allocate(tfp(npsest),cfp(npsest)) allocate(wf(npsi)) wf(1:npsi-1)=one wf(npsi)=1.0e2_wp_ - call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, & + call curfit(iopt,npsi,data%psinr,data%fpol,wf,zero,one,kspl,params%ssplf,npsest,nsf, & tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier) - call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier) -! set vacuum value used outside 0<=psin<=1 range + call splev(tfp,nsf,cfp,kspl,data%psinr(npsi:npsi),fpoli,1,ier) + ! Set vacuum value used outside 0<=data%psin<=1 range fpolas=fpoli(1) sgnbphi=sign(one,fpolas) -! free temporary arrays + ! Free temporary arrays deallocate(wrk,iwrk,wf) -! re-normalize psi after spline computation + ! Re-normalize psi after spline computation -! start with un-corrected psi - - psia=psiwbrad + ! Start with un-corrected psi + psia=data%psia psinop=0.0_wp_ psiant=1.0_wp_ -! use provided boundary to set an initial guess for the search of O/X points - - nbnd=0 - if(present(rbnd).and.present(zbnd)) then - nbnd=min(size(rbnd),size(zbnd)) - end if + ! Use provided boundary to set an initial guess + ! for the search of O/X points + nbnd=min(size(data%rbnd), size(data%zbnd)) if (nbnd>0) then - call vmaxmini(zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) - rbinf=rbnd(ibinf) - rbsup=rbnd(ibsup) - call vmaxmin(rbnd,nbnd,rbmin,rbmax) + call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) + rbinf=data%rbnd(ibinf) + rbsup=data%rbnd(ibsup) + call vmaxmin(data%rbnd,nbnd,rbmin,rbmax) else - zbinf=zv(2) - zbsup=zv(nz-1) - rbinf=rv((nr+1)/2) + zbinf=data%zv(2) + zbsup=data%zv(nz-1) + rbinf=data%rv((nr+1)/2) rbsup=rbinf - rbmin=rv(2) - rbmax=rv(nr-1) + rbmin=data%rv(2) + rbmax=data%rv(nr-1) end if -! search for exact location of the magnetic axis - - if(present(rax)) then - rax0=rax - else - rax0=0.5_wp_*(rbmin+rbmax) - end if - if(present(zax)) then - zax0=zax - else - zax0=0.5_wp_*(zbinf+zbsup) - end if + ! Search for exact location of the magnetic axis + rax0=data%rax + zax0=data%zax call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info) print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp -! search for X-point if ixp not = 0 + ! search for X-point if params%ixp /= 0 - if(present(ixp)) then - ixploc=ixp - else - ixploc=0 - end if + ixploc = params%ixp if(ixploc/=0) then if(ixploc<0) then call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) @@ -543,11 +569,13 @@ contains if (ixploc==0) then psinop=psinoptmp psiant=one-psinop -! find upper horizontal tangent point + + ! Find upper horizontal tangent point call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) zbsup=z1 rbsup=r1 -! find lower horizontal tangent point + + ! Find lower horizontal tangent point call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) zbinf=z1 rbinf=r1 @@ -555,29 +583,23 @@ contains end if print*,' ' -! save Bt value on axis (required in flux_average and used in Jcd def) -! and vacuum value B0 at ref. radius R0 (used in Jcd_astra def) + ! Save Bt value on axis (required in flux_average and used in Jcd def) + ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def) call equinum_fpol(0.0_wp_,btaxis) - btaxis=btaxis/rmaxis - if(present(r0)) then - btrcen=fpolas/r0 - rcen=r0 - else - btrcen=fpolas/rmaxis - rcen=rmaxis - end if - print'(a,f8.4)','BT_centr= ',btrcen - print'(a,f8.4)','BT_axis = ',btaxis + btaxis = btaxis/rmaxis + btrcen = fpolas/data%rvac + rcen = data%rvac + print '(a,f8.4)', 'BT_centr=', btrcen + print '(a,f8.4)', 'BT_axis =', btaxis -! compute rho_pol/rho_tor mapping based on input q profile - call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) - call set_rhospl(sqrt(psinr),rhotn) + ! 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 - subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & tx,nknt_x,ty,nknt_y,coeff,ierr) use const_and_precisions, only : wp_, comp_eps diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index bb54067..024b9c4 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -4,8 +4,8 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & 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,rtrparam_type,hcdparam_type,antctrl_type,& - eqparam_type,prfparam_type,outparam_type + 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 @@ -29,12 +29,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir real(wp_), dimension(5) :: rlim,zlim logical, save :: firstcall=.true. - type(rtrparam_type) :: rtrp - type(hcdparam_type) :: hcdp - type(antctrl_type) :: antp - type(eqparam_type) :: eqp - type(prfparam_type) :: prfp - type(outparam_type) :: outp + type(raytracing) :: rtrp + type(ecrh_cd) :: hcdp + type(antenna) :: antp + type(equilibrium) :: eqp + type(profiles) :: prfp + type(output) :: outp ! if first call read parameters from external file if (firstcall) then diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 666a7bc..c4a21f9 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -1,181 +1,246 @@ module gray_params + use const_and_precisions, only : wp_ implicit none - integer, parameter :: lenfnm=256 - integer, parameter :: headw=132,headl=21 + integer, parameter :: lenfnm = 256 + integer, parameter :: headw = 132, headl = 21 - type antctrl_type - real(wp_) :: alpha, beta, power - real(wp_) :: psi, chi - integer :: iox - integer :: ibeam - character(len=lenfnm) :: filenm - end type antctrl_type + ! Antenna/wave launcher parameters + type antenna_parameters + ! From gray_params.data: + real(wp_) :: alpha, beta ! Launching angles + real(wp_) :: power ! Initial power + real(wp_) :: psi, chi ! Initial polarisation angles + integer :: iox ! Initial wave mode + integer :: ibeam ! Beam kind + character(len=lenfnm) :: filenm ! beamdata.txt filename - type eqparam_type + ! From beamdata.txt: + real(wp_) :: fghz ! EC frequency + real(wp_) :: pos(3) ! Launcher position (tokamak frame) + real(wp_) :: w(2) ! Beam waists w(z) for the amplitude (local frame) + real(wp_) :: ri(2) ! Beam inverse radii 1/R(z) for the phase (local frame) + real(wp_) :: phi(2) ! Axes orientation angles for amplitude, phase (local frame) + end type + + ! MHD equilibrium parameters + type equilibrium_parameters real(wp_) :: ssplps, ssplf, factb integer :: sgnb, sgni, ixp integer :: iequil, icocos, ipsinorm, idesc, ifreefmt character(len=lenfnm) :: filenm - end type eqparam_type + end type - type prfparam_type - real(wp_) :: psnbnd, sspld, factne, factte - integer :: iscal, irho !, icrho, icte, icne, iczf + ! Kinetic plasma profiles + type profiles_parameters + real(wp_) :: psnbnd ! plasma boundary (ψ: ne(ψ)=0) + real(wp_) :: sspld, factne, factte + integer :: iscal, irho integer :: iprof character(len=lenfnm) :: filenm - end type prfparam_type + end type - type rtrparam_type + ! Raytracing parameters + type raytracing_parameters real(wp_) :: rwmax, dst integer :: nrayr, nrayth, nstep integer :: igrad, idst, ipass, ipol - end type rtrparam_type + end type - type hcdparam_type + ! EC resonant heating & Current Drive parameters + type ecrh_cd_parameters integer :: iwarm, ilarm, imx, ieccd - end type hcdparam_type + end type - type outparam_type - integer :: ipec, nrho, istpr, istpl - end type outparam_type + ! Output data parameters + type output_parameters + integer :: ipec ! Grid spacing for profiles (profili EC) + integer :: nrho ! Number of grid steps for EC profiles + 1 + integer :: istpr ! Subsampling factor for beam cross-section (units 8, 12) + integer :: istpl ! Subsampling factor for outer rays (unit 33) + end type - integer, save :: iequil,iprof,ipol - integer, save :: iwarm,ilarm,imx,ieccd - integer, save :: igrad,idst,ipass - integer, save :: istpr0,istpl0 - integer, save :: ipec,nnd + ! Other parameters + type misc_parameters + real(wp_) :: rwall ! R of the limiter (fallback) + end type + + ! MHD equilibrium data + type equilibrium_data + real(wp_), allocatable :: rv(:) ! R of the uniform grid + real(wp_), allocatable :: zv(:) ! Z of the uniform grid + real(wp_), allocatable :: rlim(:) ! R of the limiter contour (wall) + real(wp_), allocatable :: zlim(:) ! Z of the limiter contour + real(wp_), allocatable :: rbnd(:) ! R of the boundary contour (plasma) + real(wp_), allocatable :: zbnd(:) ! Z of the boundary contour + real(wp_), allocatable :: fpol(:) ! Poloidal current function + real(wp_), allocatable :: qpsi(:) ! Safety factor on the flux grid + real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid + real(wp_), allocatable :: psinr(:) ! Poloidal flux + real(wp_) :: psia ! Poloidal flux at edge - flux at magnetic axis + real(wp_) :: rvac ! Reference R₀ (B = B₀R₀ without the plasma) + real(wp_) :: rax ! R of the magnetic axis + real(wp_) :: zax ! Z of the magnetic axis + end type + + ! Kinetic plasma profiles data + type profiles_data + real(wp_), allocatable :: psrad(:) ! Radial coordinate + real(wp_), allocatable :: terad(:) ! Electron temperature profile + real(wp_), allocatable :: derad(:) ! Electron density profile + real(wp_), allocatable :: zfc(:) ! Effective charge profile + end type + + ! All GRAY parameters + type gray_parameters + type(antenna_parameters) :: antenna + type(equilibrium_parameters) :: equilibrium + type(profiles_parameters) :: profiles + type(raytracing_parameters) :: raytracing + type(ecrh_cd_parameters) :: ecrh_cd + type(output_parameters) :: output + type(misc_parameters) :: misc + end type + + ! All GRAY input data + type gray_data + type(equilibrium_data) :: equilibrium + type(profiles_data) :: profiles + end type + + ! GRAY final results (only some) + type gray_results + real(wp_) :: pabs ! Total absorbed power + real(wp_) :: icd ! Total driven current + real(wp_), allocatable :: dpdv(:) ! Absorbed power density + real(wp_), allocatable :: jcd(:) ! Driven current density + end type + + ! Global variables exposed for graycore + integer, save :: iequil, iprof, ipol + integer, save :: iwarm, ilarm, imx, ieccd + integer, save :: igrad, idst, ipass + integer, save :: istpr0, istpl0 + integer, save :: ipec, nnd contains - subroutine print_params(rtrparam,hcdparam,antctrl,eqparam,rwall, & - prfparam,outparam,strout) + subroutine print_parameters(params, strout) implicit none - ! arguments - type(rtrparam_type), intent(in) :: rtrparam - type(hcdparam_type), intent(in) :: hcdparam - type(antctrl_type), intent(in) :: antctrl - type(eqparam_type), intent(in) :: eqparam - real(wp_), intent(in) :: rwall - type(prfparam_type), intent(in) :: prfparam - type(outparam_type), intent(in) :: outparam - character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21) + ! subroutine arguments + type(gray_parameters), intent(in) :: params + character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21) ! local variables - character(len=8) :: rdat - character(len=10) :: rtim + character(len=8) :: date + character(len=10) :: time #ifndef REVISION character(len=*), parameter :: REVISION="unknown" #endif - ! date and time - call date_and_time(rdat,rtim) - write(strout(1),'("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') & - rdat(1:4),rdat(5:6),rdat(7:8),rtim(1:2),rtim(3:4),rtim(5:10) + ! Date and time + call date_and_time(date, time) + write(strout(1), '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') & + date(1:4), date(5:6), date(7:8), & + time(1:2), time(3:4), time(5:10) ! Git revision - write(strout(2),'("# GRAY Git revision: ",a)') REVISION + write(strout(2), '("# GRAY Git revision: ",a)') REVISION - ! equilibrium input data - if (eqparam%iequil > 0) then - write(strout(3),'("# EQL input: ",a)') trim(eqparam%filenm) - !!! missing values - write(strout(7),'("# EQL B0 R0 aminor Rax zax:",5(1x,e12.5))') & + ! MHD equilibrium parameters + if (params%equilibrium%iequil > 0) then + write(strout(3), '("# EQL input: ",a)') trim(params%equilibrium%filenm) + ! TODO add missing values + write(strout(7), '("# EQL B0 R0 aminor Rax zax:",5(1x,e12.5))') & 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_ else - write(strout(3),'("# EQL input: N/A (vacuum)")') - write(strout(7),'("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') + write(strout(3), '("# EQL input: N/A (vacuum)")') + write(strout(7), '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') end if - write(strout(4),'("# EQL iequil sgnb sgni factb:",3(1x,i4),1x,e12.5)') & - eqparam%iequil, eqparam%sgnb, eqparam%sgni, eqparam%factb - if (eqparam%iequil > 1) then - write(strout(5),'("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,i4))') & - eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt - write(strout(6),'("# EQL ssplps ssplf ixp:",2(1x,e12.5),1x,i4)') & - eqparam%ssplps, eqparam%ssplf, eqparam%ixp + write(strout(4), '("# EQL iequil sgnb sgni factb:",3(1x,i4),1x,e12.5)') & + params%equilibrium%iequil, params%equilibrium%sgnb, params%equilibrium%sgni, params%equilibrium%factb + if (params%equilibrium%iequil > 1) then + write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,i4))') & + params%equilibrium%icocos, params%equilibrium%ipsinorm, params%equilibrium%idesc, params%equilibrium%ifreefmt + write(strout(6), '("# EQL ssplps ssplf ixp:",2(1x,e12.5),1x,i4)') & + params%equilibrium%ssplps, params%equilibrium%ssplf, params%equilibrium%ixp else - write(strout(5),'("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")') - write(strout(6),'("# EQL ssplps ssplf ixp: N/A (analytical)")') + write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")') + write(strout(6), '("# EQL ssplps ssplf ixp: N/A (analytical)")') end if - ! profiles input data - if (eqparam%iequil > 0) then - write(strout(8),'("# PRF input: ",a)') trim(prfparam%filenm) - write(strout(9),'("# PRF iprof iscal factne factte:",2(1x,i4),2(1x,e12.5))') & - prfparam%iprof,prfparam%iscal,prfparam%factne,prfparam%factte - if (prfparam%iprof > 0) then - write(strout(10),'("# PRF irho psnbnd sspld:",1x,i4,2(1x,e12.5))') & - prfparam%irho,prfparam%psnbnd,prfparam%sspld + ! Profiles parameters + if (params%equilibrium%iequil > 0) then + write(strout(8), '("# PRF input: ",a)') trim(params%profiles%filenm) + write(strout(9), '("# PRF iprof iscal factne factte:",2(1x,i4),2(1x,e12.5))') & + params%profiles%iprof, params%profiles%iscal, params%profiles%factne, params%profiles%factte + if (params%profiles%iprof > 0) then + write(strout(10), '("# PRF irho psnbnd sspld:",1x,i4,2(1x,e12.5))') & + params%profiles%irho,params%profiles%psnbnd,params%profiles%sspld else - write(strout(10),'("# PRF irho psnbnd sspld: N/A (analytical)")') + write(strout(10), '("# PRF irho psnbnd sspld: N/A (analytical)")') end if - !!! missing values - write(strout(11),'("# PRF Te0 ne0 Zeff0:",3(1x,e12.5))') & + ! TODO: add missing values + write(strout(11), '("# PRF Te0 ne0 Zeff0:",3(1x,e12.5))') & 0._wp_, 0._wp_, 0._wp_ else - write(strout(8),'("# PRF input: N/A (vacuum)")') - write(strout(9),'("# PRF iprof iscal factne factte: N/A (vacuum)")') - write(strout(10),'("# PRF irho psnbnd sspld: N/A (vacuum)")') - write(strout(11),'("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') + write(strout(8), '("# PRF input: N/A (vacuum)")') + write(strout(9), '("# PRF iprof iscal factne factte: N/A (vacuum)")') + write(strout(10), '("# PRF irho psnbnd sspld: N/A (vacuum)")') + write(strout(11), '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') end if - ! launch parameters - write(strout(12),'("# ANT input: ",a)') trim(antctrl%filenm) - write(strout(13),'("# ANT ibeam iox psi chi:",2(1x,i4),2(1x,e12.5))') & - antctrl%ibeam, antctrl%iox, antctrl%psi, antctrl%chi - write(strout(14),'("# ANT alpha beta power:",3(1x,e12.5))') & - antctrl%alpha, antctrl%beta, antctrl%power - !!! missing values - write(strout(15),'("# ANT x0 y0 z0:",3(1x,e12.5))') & + ! Antenna parameters + write(strout(12), '("# ANT input: ",a)') trim(params%antenna%filenm) + write(strout(13), '("# ANT ibeam iox psi chi:",2(1x,i4),2(1x,e12.5))') & + params%antenna%ibeam, params%antenna%iox, params%antenna%psi, params%antenna%chi + write(strout(14), '("# ANT alpha beta power:",3(1x,e12.5))') & + params%antenna%alpha, params%antenna%beta, params%antenna%power + ! TODO: add missing values + write(strout(15), '("# ANT x0 y0 z0:",3(1x,e12.5))') & 0._wp_, 0._wp_, 0._wp_ - !!! missing values - write(strout(16),'("# ANT wx wy Rcix Rciy psiw psir:",6(1x,e12.5))') & + ! TODO: add missing values + write(strout(16), '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,e12.5))') & 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_ - ! wall parameters - write(strout(17),'("# RFL rwall:",1x,e12.5)') rwall + ! Other parameters + write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall ! code parameters - write(strout(18),'("# COD igrad idst ipass ipol:",4(1x,i4))') & - rtrparam%igrad, rtrparam%idst, rtrparam%ipass, rtrparam%ipol - write(strout(19),'("# COD nrayr nrayth nstep rwmax dst:",3(1x,i4),2(1x,e12.5))') & - rtrparam%nrayr, rtrparam%nrayth, rtrparam%nstep, rtrparam%rwmax, rtrparam%dst - write(strout(20),'("# COD iwarm ilarm imx ieccd:",4(1x,i4))') & - hcdparam%iwarm, hcdparam%ilarm, hcdparam%imx, hcdparam%ieccd - write(strout(21),'("# COD ipec nrho istpr istpl:",4(1x,i4))') & - outparam%ipec, outparam%nrho, outparam%istpr, outparam%istpl - end subroutine print_params + write(strout(18), '("# COD igrad idst ipass ipol:",4(1x,i4))') & + params%raytracing%igrad, params%raytracing%idst, params%raytracing%ipass, params%raytracing%ipol + write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",3(1x,i4),2(1x,e12.5))') & + params%raytracing%nrayr, params%raytracing%nrayth, params%raytracing%nstep, params%raytracing%rwmax, params%raytracing%dst + write(strout(20), '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') & + params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx, params%ecrh_cd%ieccd + write(strout(21), '("# COD ipec nrho istpr istpl:",4(1x,i4))') & + params%output%ipec, params%output%nrho, params%output%istpr, params%output%istpl + end subroutine print_parameters - subroutine read_params(filenm,rtrparam,hcdparam,antctrl,eqparam,rwall, & - prfparam,outparam,unit) + subroutine read_parameters(filename, params, unit) use utils, only : get_free_unit implicit none - ! arguments - character(len=*), intent(in) :: filenm - type(rtrparam_type), intent(out) :: rtrparam - type(hcdparam_type), intent(out) :: hcdparam - type(antctrl_type), intent(out) :: antctrl - type(eqparam_type), intent(out) :: eqparam - real(wp_), intent(out) :: rwall - type(prfparam_type), intent(out) :: prfparam - type(outparam_type), intent(out) :: outparam + ! subrouting arguments + character(len=*), intent(in) :: filename + type(gray_parameters), intent(out) :: params integer, intent(in), optional :: unit ! local variables integer :: u, iostat if (present(unit)) then - u=unit + u = unit else u = get_free_unit() end if - open(u,file=filenm,status='old',action='read',iostat=iostat) - if (iostat>0) then - print'(a)', 'gray_params file not found!' + + open(u, file=filename, status='old', action='read', iostat=iostat) + if (iostat > 0) then + print '(a)', 'gray_params file not found!' call EXIT(1) end if @@ -185,18 +250,18 @@ contains ! nrayth :number of rays in angular direction ! rwmax :normalized maximum radius of beam power ! rwmax=1 -> :last ray at radius = waist - read(u,*) rtrparam%nrayr, rtrparam%nrayth, rtrparam%rwmax + read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, params%raytracing%rwmax ! igrad=0 :optical ray-tracing, initial conditions as for beam ! igrad=1 :quasi-optical ray-tracing ! igrad=-1 :ray-tracing, init. condit. ! from center of mirror and with angular spread ! ipass=1/2 :1 or 2 passes into plasma ! ipol=0 :compute mode polarization at antenna, ipol=1 use polariz angles - read(u,*) rtrparam%igrad, rtrparam%ipass, rtrparam%ipol + read(u, *) params%raytracing%igrad, params%raytracing%ipass, params%raytracing%ipol ! dst :integration step ! nstep :maximum number of integration steps ! idst=0/1/2 :0 integration in s, 1 integr. in ct, 2 integr. in Sr - read(u,*) rtrparam%dst, rtrparam%nstep, rtrparam%idst + read(u, *) params%raytracing%dst, params%raytracing%nstep, params%raytracing%idst ! Heating & Current drive ! ======================================================================== @@ -207,104 +272,102 @@ contains ! ilarm :order of larmor expansion ! imx :max n of iterations in dispersion, imx<0 uses 1st ! iteration in case of failure after |imx| iterations - read(u,*) hcdparam%iwarm,hcdparam%ilarm,hcdparam%imx + read(u, *) params%ecrh_cd%iwarm,params%ecrh_cd%ilarm,params%ecrh_cd%imx ! ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models - read(u,*) hcdparam%ieccd + read(u, *) params%ecrh_cd%ieccd ! Wave launcher ! ======================================================================== ! alpha0, beta0 (cartesian) launching angles - read(u,*) antctrl%alpha, antctrl%beta + read(u, *) params%antenna%alpha, params%antenna%beta ! p0mw injected power (MW) - read(u,*) antctrl%power + read(u, *) params%antenna%power ! abs(iox)=1/2 OM/XM ! psipol0,chipol0 polarization angles at the antenna (if iox<0) - read(u,*) antctrl%iox, antctrl%psi, antctrl%chi - ! ibeam=0 :read data for beam as above - ! ibeam=1 :read data from file simple astigmatic beam - ! ibeam=2 :read data from file general astigmatic beam - read(u,*) antctrl%ibeam - read(u,*) antctrl%filenm + read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi + read(u, *) params%antenna%ibeam + read(u, *) params%antenna%filenm - ! Magnetic equilibrium + ! MHD equilibrium ! ======================================================================== - ! iequil=0 :vacuum + ! iequil=0 :vacuum (no plasma) ! iequil=1 :analytical equilibrium ! iequil=2 :read eqdsk - read(u,*) eqparam%iequil - read(u,*) eqparam%filenm + ! iequil=2 :read eqdsk, data only valid inside last closed flux surface + read(u, *) params%equilibrium%iequil + read(u, *) params%equilibrium%filenm ! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012 ! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004 - read(u,*) eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt + read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, params%equilibrium%idesc, params%equilibrium%ifreefmt ! ixp=0,-1,+1 : no X point , bottom/up X point ! ssplps : spline parameter for psi interpolation - read(u,*) eqparam%ixp, eqparam%ssplps !, eqparam%ssplf - eqparam%ssplf=0.01_wp_ + read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps !, params%equilibrium%ssplf + params%equilibrium%ssplf=0.01_wp_ ! signum of toroidal B and I ! factb factor for magnetic field (only for numerical equil) ! scaling adopted: beta=const, qpsi=const, nustar=const - read(u,*) eqparam%sgnb, eqparam%sgni, eqparam%factb + read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, params%equilibrium%factb ! Wall ! ======================================================================== - read(u,*) rwall + read(u, *) params%misc%rwall ! Profiles ! ======================================================================== ! iprof=0 :analytical density and temp. profiles ! iprof>0 :numerical density and temp. profiles - read(u,*) prfparam%iprof, prfparam%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin - read(u,*) prfparam%filenm + read(u, *) params%profiles%iprof, params%profiles%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin + read(u, *) params%profiles%filenm ! psbnd value of psi ( > 1 ) of density boundary - read(u,*) prfparam%psnbnd, prfparam%sspld + read(u, *) params%profiles%psnbnd, params%profiles%sspld ! prfparam%sspld=0.001_wp_ ! iscal :ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling ! factT factn :factor for Te&ne scaling - read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal + read(u, *) params%profiles%factte, params%profiles%factne, params%profiles%iscal ! Output ! ======================================================================== ! ipec=0/1 :pec profiles grid in psi/rhop ! nrho :number of grid steps for pec profiles +1 - read(u,*) outparam%ipec, outparam%nrho + read(u, *) params%output%ipec, params%output%nrho ! istpr0 :projection step = dsdt*istprj ! istpl0 :plot step = dsdt*istpl - read(u,*) outparam%istpr, outparam%istpl + read(u, *) params%output%istpr, params%output%istpl + close(u) - end subroutine read_params + end subroutine read_parameters - subroutine set_codepar(eqparam,prfparam,outparam,rtrparam,hcdparam) + subroutine set_globals(params) implicit none - type(eqparam_type), intent(in) :: eqparam - type(prfparam_type), intent(in) :: prfparam - type(outparam_type), intent(in) :: outparam - type(rtrparam_type), intent(in) :: rtrparam - type(hcdparam_type), intent(in) :: hcdparam - iequil=eqparam%iequil - iprof=prfparam%iprof + ! subroutine arguments + type(gray_parameters), intent(in) :: params - ipec=outparam%ipec - nnd=outparam%nrho - istpr0=outparam%istpr - istpl0=outparam%istpl + iequil = params%equilibrium%iequil + iprof = params%profiles%iprof - ipol=rtrparam%ipol - igrad=rtrparam%igrad - idst=rtrparam%idst - ipass=rtrparam%ipass - if (rtrparam%nrayr<5) then - igrad=0 - print*,' nrayr < 5 ! => OPTICAL CASE ONLY' - print*,' ' + ipec = params%output%ipec + nnd = params%output%nrho + istpr0 = params%output%istpr + istpl0 = params%output%istpl + + ipol = params%raytracing%ipol + igrad = params%raytracing%igrad + idst = params%raytracing%idst + ipass = params%raytracing%ipass + + if (params%raytracing%nrayr < 5) then + igrad = 0 + print *, ' nrayr < 5 ! => OPTICAL CASE ONLY' + print *, ' ' end if - iwarm=hcdparam%iwarm - ilarm=hcdparam%ilarm - imx=hcdparam%imx - ieccd=hcdparam%ieccd + iwarm = params%ecrh_cd%iwarm + ilarm = params%ecrh_cd%ilarm + imx = params%ecrh_cd%imx + ieccd = params%ecrh_cd%ieccd - end subroutine set_codepar + end subroutine set_globals end module gray_params diff --git a/src/graycore.f90 b/src/graycore.f90 index 77a9d5b..fa955a6 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -4,56 +4,41 @@ module graycore contains - subroutine gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, & - eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & - p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & - psipol0,chipol0, dpdv,jcd,pabs,icd, outp,rtrp,hcdp,ierr, rhout) + subroutine gray_main(params, data, results, error, rhout) use const_and_precisions, only : zero, one, degree, comp_tiny use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl - use dispersion, only : expinit - use gray_params, only : eqparam_type, prfparam_type, outparam_type, & - rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, & - iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl, ipass - use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff - use beamdata, only : pweight, rayi2jk - use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & - zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q, psnbd - use errcodes, only : check_err, print_errn, print_errhcd + use dispersion, only : expinit + use gray_params, only : gray_parameters, gray_data, gray_results, print_parameters, & + iwarm, ipec, istpr0, igrad, headw, headl, ipass + use beams, only : xgygcoeff, launchangles2n + use beamdata, only : pweight, rayi2jk + use equilibrium, only : unset_eqspl, unset_rhospl, unset_q + use errcodes, only : check_err, print_errn, print_errhcd use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst - use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & - rhop_tab, rhot_tab - use limiter, only : set_lim, unset_lim - use utils, only : vmaxmin - use reflections, only : inside - use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & - initmultipass, turnoffray, plasma_in, plasma_out, wall_out - use units, only : ucenr + use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst + use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & + rhop_tab, rhot_tab + use limiter, only : limiter_unset_globals=>unset_globals + use utils, only : vmaxmin + use reflections, only : inside + use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & + initmultipass, turnoffray, plasma_in, plasma_out, wall_out + use units, only : ucenr implicit none -! arguments - type(eqparam_type), intent(in) :: eqp - type(prfparam_type), intent(in) :: prfp - type(outparam_type), intent(in) :: outp - type(rtrparam_type), intent(in) :: rtrp - type(hcdparam_type), intent(in) :: hcdp - real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc - real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi - real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim - real(wp_), dimension(:,:), intent(in) :: psin - real(wp_), intent(in) :: psia, rvac, rax, zax - integer, intent(in) :: iox0 - real(wp_), intent(in) :: p0, fghz, psipol0, chipol0 - real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir - real(wp_), dimension(3), intent(in) :: xv0 + ! Subroutine arguments + type(gray_parameters), intent(in) :: params + type(gray_data), intent(in) :: data + type(gray_results), intent(out) :: results - real(wp_), intent(out) :: pabs,icd - real(wp_), dimension(:), intent(out) :: dpdv,jcd + ! Predefined grid for the output profiles (optional) real(wp_), dimension(:), intent(in), optional :: rhout - integer, intent(out) :: ierr -! local variables + ! Exit code + integer, intent(out) :: error + + ! local variables real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) character, dimension(2), parameter :: mode=(/'O','X'/) @@ -66,11 +51,17 @@ contains real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0 real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg - real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() + + ! Ray variables + real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() - integer :: i,j,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt + + ! i: integration step, jk: global ray index + integer :: i, jk + + integer :: iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt integer :: ip,ib,iopmin,ipar,iO - integer :: igrad_b,ipol,istop_pass,nbeam_pass,nlim + integer :: igrad_b,istop_pass,nbeam_pass,nlim logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null() @@ -88,101 +79,92 @@ contains logical, dimension(:), pointer :: iwait=>null() logical, dimension(:,:), pointer :: iroff=>null() -! parameters log in file headers + ! parameters log in file headers character(len=headw), dimension(headl) :: strheader - type(antctrl_type) :: antp - real(wp_) :: rwall -! ======== set environment BEGIN ======== - call set_codepar(eqp,prfp,outp,rtrp,hcdp) + ! ======== set environment BEGIN ======== + ! Number of limiter contourn points + nlim = size(data%equilibrium%zlim) - call set_lim(rlim,zlim) - nlim = size(zlim) + ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) + call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) - if(iequil<2) then - call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) - else - call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, & - rvac, rax,zax, rbnd,zbnd, eqp%ixp) - ! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol) - end if -! compute flux surface averaged quantities - call flux_average ! requires frhotor for dadrhot,dvdrhot + ! Compute the initial cartesian wavevector (anv0) + ! from launch angles α,β and the position x₀: + ! NR(α, β, x₀) + ! Nφ(α, β, x₀) + ! Nz(α, β, x₀) + call launchangles2n(params%antenna, anv0) - if(iprof==0) then - call set_prfan(terad,derad,zfc) - else - call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd) - end if - - call xgygcoeff(fghz,ak0,bres,xgcn) - call launchangles2n(alpha0,beta0,xv0,anv0) - call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) + ! Initialise the ray variables (beamtracing) + call init_btr(params%raytracing, yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, & + tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv) + ! Initialise the dispersion module if(iwarm > 1) call expinit + + ! Initialise the magsurf_data module + call flux_average ! requires frhotor for dadrhot,dvdrhot - call pec_init(ipec,rhout) - nnd=size(rhop_tab) - call alloc_multipass(nnd,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & - stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) + ! Initialise the output profiles + call pec_init(ipec, rhout) + nnd = size(rhop_tab) ! number of radial profile points -! ========= set environment END ========= + call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & + stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & + pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) -! ======== pre-proc prints BEGIN ======== - antp%alpha=alpha0 - antp%beta=beta0 - antp%power=p0 - antp%psi=psipol0 - antp%chi=chipol0 - antp%iox=iox0 -!!!!! missing values - antp%ibeam=0 - antp%filenm='' - rwall=0._wp_ - psnbd=prfp%psnbnd ! plasma boundary - ipol=rtrp%ipol + ! Allocate memory for the results... + allocate(results%dpdv(params%output%nrho)) + allocate(results%jcd(params%output%nrho)) - pabs=zero ! gray_main output - icd=zero - dpdv=zero - jcd=zero - - call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) - call print_headers(strheader,0) -! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout - call print_surfq((/1.5_wp_,2.0_wp_/)) -! print - print*,' ' - print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 - print'(a,4f8.3)','x00, y00, z00 = ',xv0 -! print Btot=Bres -! print ne, Te, q, Jphi versus psi, rhop, rhot + ! ...and initialise them + results%pabs = zero + results%icd = zero + results%dpdv = zero + results%jcd = zero + ! ========= set environment END ========= + + ! ======== pre-proc prints BEGIN ======== + call print_parameters(params, strheader) + call print_headers(strheader) + + ! print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout + call print_surfq([1.5_wp_, 2.0_wp_]) + + ! print initial position + print *, '' + print '(a,2f8.3)', 'alpha0, beta0 = ', params%antenna%alpha, params%antenna%beta + print '(a,4f8.3)', 'x00, y00, z00 = ', params%antenna%pos + + ! print Btot=Bres + ! print ne, Te, q, Jphi versus psi, rhop, rhot call print_bres(bres) call print_prof - call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2),sin(beta0*degree)) + call print_maps(bres, xgcn, & + 0.01_wp_*sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2), & + sin(params%antenna%beta*degree)) + ! ========= pre-proc prints END ========= -! ========= pre-proc prints END ========= - -! =========== main loop BEGIN =========== - call initmultipass(ipol,iox0,iroff,yynext,yypnext,yw0,ypw0, & - stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) + ! =========== main loop BEGIN =========== + call initmultipass(params%raytracing%ipol, params%antenna%iox, & + iroff,yynext,yypnext,yw0,ypw0, & + stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) - if(ipol.eq.0) then - if(iox0.eq.2) then ! only X mode on 1st pass + if(params%raytracing%ipol .eq. 0) then + if(params%antenna%iox .eq. 2) then ! only X mode on 1st pass cpl0 = (/zero,one/) else ! only O mode on 1st pass cpl0 = (/one,zero/) end if end if - + sox=one ! mode inverted for each beam iox=2 ! start with O: sox=-1, iox=1 - psipol=psipol0 - chipol=chipol0 - call pweight(p0,p0jk) + psipol = params%antenna%psi + chipol = params%antenna%chi + call pweight(params%antenna%power, p0jk) nbeam_pass=1 ! max n of beam per pass index_rt=0 ! global beam index: 1,O 2,X 1st pass @@ -219,7 +201,6 @@ contains end if call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) -! call print_headers((/' '/),index_rt) if(ip.eq.1) then ! 1st pass igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass @@ -230,14 +211,19 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions + call ic_gb(params%antenna%pos, anv0, ak0, & + params%antenna%w(1),params%antenna%w(2), & + params%antenna%ri(1),params%antenna%ri(2), & + params%antenna%phi(1),params%antenna%phi(2), & + yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization do jk=1,nray zzm = yw(3,jk)*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ - if(inside(rlim,zlim,nlim,rrm,zzm)) then ! * start propagation in/outside vessel? + if(inside(data%equilibrium%rlim, data%equilibrium%zlim, & + nlim, rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside else iow(jk) = 0 ! + outside @@ -273,7 +259,7 @@ contains ! update position and grad if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) - ierr = 0 + error = 0 istop = 0 ! stop flag for current beam iopmin = 10 @@ -289,7 +275,7 @@ contains ierrn,igrad_b) ! update global error code and print message if(ierrn/=0) then - ierr = ior(ierr,ierrn) + error = ior(error,ierrn) call print_errn(ierrn,i,anpl) end if @@ -297,8 +283,9 @@ contains zzm = xv(3)*0.01_wp_ rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ - ins_pl = (psinv>=zero .and. psinv=zbinf .and. zzm<=zbsup ! in/out plasma? - ins_wl = inside(rlim,zlim,nlim,rrm,zzm) ! in/out vessel? + ins_pl = (psinv>=zero .and. psinv continue current pass - if(ipol.eq.0) then ! + IF single mode propagation + if(params%raytracing%ipol .eq. 0) then ! + IF single mode propagation cpl = cpl0 p0ray(jk) = p0ray(jk)*cpl(iox) else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays @@ -376,7 +363,7 @@ contains psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, & ierrn,igrad_b) ! * update derivatives after reflection if(ierrn/=0) then ! * update global error code and print message - ierr = ior(ierr,ierrn) + error = ior(error,ierrn) call print_errn(ierrn,i,anpl) end if @@ -429,7 +416,7 @@ contains call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) if(ierrhcd/=0) then - ierr = ior(ierr,ierrhcd) + error = ior(error,ierrhcd) call print_errhcd(ierrhcd,i,anprre,anprim,alpha) end if else @@ -464,9 +451,6 @@ contains call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, & btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) -! call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, & -! btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & -! nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) end if end do @@ -484,7 +468,7 @@ contains end if ! check for any error code and stop if necessary - call check_err(ierr,istop) + call check_err(error,istop) ! test whether further trajectory integration is unnecessary call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling ! if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam @@ -563,8 +547,8 @@ contains ! ============ beam loop END ============ ! ======= cumulative prints BEGIN ======= - pabs = pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] - icd = icd + sum(icd_pass) + results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] + results%icd = results%icd + sum(icd_pass) ! print final results for pass on screen write(*,*) @@ -581,16 +565,11 @@ contains ! print final results on screen write(*,*) - write(*,'(a)') '## Final results:' - write(*,'(a,f9.4)') '## Pabs_tot (MW) = ',pabs - write(*,'(a,f9.4)') '## I_tot (kA) = ',icd*1.0e3_wp_ - + write(*,'(a)') '## Final results:' + write(*,'(a,f9.4)') '## Pabs_tot (MW) = ', results%pabs + write(*,'(a,f9.4)') '## I_tot (kA) = ', results%icd*1.0e3_wp_ + ! ========== free memory BEGIN ========== - call unset_eqspl - call unset_q - call unset_rhospl - call unset_prfspl - call unset_lim call dealloc_surfvec call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) @@ -603,170 +582,6 @@ contains - subroutine sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, & - eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & - p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & - psipol0,chipol0, jphi,jcd,dpdv,currins,pins,pabs,icd, & - jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & - outp,rtrp,hcdp,ierr, rhout) - use const_and_precisions, only : zero, one, degree - use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl - use dispersion, only : expinit - use gray_params, only : eqparam_type, prfparam_type, outparam_type, & - rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, & - iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl - use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff - use beamdata, only : pweight, rayi2jk - use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & - zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q - use errcodes, only : check_err, print_errn, print_errhcd - use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst - use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & - rhop_tab, rhot_tab - use limiter, only : set_lim, unset_lim - use utils, only : vmaxmin - implicit none -! arguments - type(eqparam_type), intent(in) :: eqp - type(prfparam_type), intent(in) :: prfp - type(outparam_type), intent(in) :: outp - type(rtrparam_type), intent(in) :: rtrp - type(hcdparam_type), intent(in) :: hcdp - - real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc - real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi - real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim - real(wp_), dimension(:,:), intent(in) :: psin - real(wp_), intent(in) :: psia, rvac, rax, zax - integer, intent(in) :: iox0 - real(wp_), intent(in) :: p0, fghz, psipol0, chipol0 - real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir - real(wp_), dimension(3), intent(in) :: xv0 - - real(wp_), intent(in) :: pabs,icd - real(wp_), dimension(:), intent(in) :: jphi,jcd,dpdv,currins,pins - real(wp_), intent(out) :: jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - - real(wp_), dimension(:), intent(in), optional :: rhout - - integer, intent(out) :: ierr -! local variables - real(wp_), parameter :: taucr = 12._wp_ - - real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre - real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0 - real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx - real(wp_) :: drhotp,drhotj,dpdvmx,jphimx - - real(wp_), dimension(3) :: xv,anv0,anv,bv - real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() - real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() - integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1 - logical :: ins_pl, somein, allout - - real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null() - real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), & - ccci0=>null() - real(wp_), dimension(:), pointer :: p0jk=>null() - complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() - integer, dimension(:), pointer :: iiv=>null() - -! parameters log in file headers - character(len=headw), dimension(headl) :: strheader - type(antctrl_type) :: antp - real(wp_) :: rwall - -! ======== set environment BEGIN ======== - call set_codepar(eqp,prfp,outp,rtrp,hcdp) - - call set_lim(rlim,zlim) - - if(iequil<2) then - call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) - else - call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, & - rvac, rax,zax, rbnd,zbnd, eqp%ixp) - ! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol) - end if -! compute flux surface averaged quantities - call flux_average ! requires frhotor for dadrhot,dvdrhot - - if(iprof==0) then - call set_prfan(terad,derad,zfc) - else - call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd) - end if - - call xgygcoeff(fghz,ak0,bres,xgcn) - call launchangles2n(alpha0,beta0,xv0,anv0) - call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - - if(iwarm > 1) call expinit - -! ======= set environment END ====== - -! ======= pre-proc prints BEGIN ====== - antp%alpha=alpha0 - antp%beta=beta0 - antp%power=p0 - antp%psi=psipol0 - antp%chi=chipol0 - antp%iox=iox0 -!!!!! missing values - antp%ibeam=0 - antp%filenm='' - rwall=0._wp_ - call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) - call print_headers(strheader, 0) -! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout - call print_surfq((/1.5_wp_,2.0_wp_/)) -! print - print*,' ' - print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 - print'(a,4f8.3)','x00, y00, z00 = ',xv0 -! print Btot=Bres -! print ne, Te, q, Jphi versus psi, rhop, rhot - call print_bres(bres) - call print_prof - call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), & - sin(beta0*degree)) -! ======= pre-proc prints END ====== - -! ======= post-proc BEGIN ====== - -! compute power and current density profiles for all rays - call pec_init(ipec,rhout) - nnd=size(rhop_tab) -! print power and current density profiles - call print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt) -! compute profiles width - call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, & - rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & - rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) -! print 0D results - call print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, & - drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, & - st,psipol,chipol,index_rt,p0,zero,zero) ! cpl1=cpl2=0 for a single pass - -! ======= post-proc END ====== - -! ======= free memory BEGIN ====== - call unset_eqspl - call unset_q - call unset_rhospl - call unset_prfspl - call unset_lim - call dealloc_surfvec - call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - call dealloc_pec -! ======= free memory END ====== - end subroutine sum_profiles - subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) use const_and_precisions, only : wp_, zero @@ -1167,7 +982,7 @@ contains subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & - bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr,igrad) + bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use errcodes, only : pnpl @@ -1181,7 +996,7 @@ contains real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm real(wp_), dimension(3), intent(out) :: bv - integer, intent(out) :: ierr + integer, intent(out) :: error real(wp_), dimension(3), intent(out) :: derxg integer, intent(in) :: igrad ! local variables @@ -1196,12 +1011,12 @@ contains call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) - ierr=0 + error=0 if( abs(anpl) > anplth1) then if(abs(anpl) > anplth2) then - ierr=ibset(ierr,pnpl+1) + error=ibset(error,pnpl+1) else - ierr=ibset(ierr,pnpl) + error=ibset(error,pnpl) end if end if end subroutine ywppla_upd @@ -1692,7 +1507,7 @@ contains subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & - sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) + sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error) use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ use gray_params, only : iwarm,ilarm,ieccd,imx use coreprofiles, only : fzeff @@ -1707,7 +1522,7 @@ contains real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox real(wp_),intent(out) :: anprre,anprim,alpha,didp integer, intent(out) :: nhmin,nhmax,iokhawa - integer, intent(out) :: ierr + integer, intent(out) :: error ! local constants real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ ! local variables @@ -1726,7 +1541,7 @@ contains nhmin=0 nhmax=0 iokhawa=0 - ierr=0 + error=0 if(tekev>zero) then ! absorption computation @@ -1734,13 +1549,13 @@ contains call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) if(nhmin.gt.0) then lrm=max(ilarm,nhmax) - call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, & + call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,error,anprre,anprim, & iwarm,imx,ex,ey,ez) akim=ak0*anprim ratiovgr=2.0_wp_*anpr/derdnm!*vgm alpha=2.0_wp_*akim*ratiovgr if(alpha= zero .and. psinv < psdbnd) .and. & -! (psinv >= one .or. (zzm >= zbinf .and. zzm <= zbsup)) -! end function inside_plasma -! -! -! -! subroutine vacuum_rt(xv0,anv0,xvend,dstvac,ivac) -! use const_and_precisions, only : wp_ -! use beamdata, only : dst -! use limiter, only : rlim,zlim,nlim -! implicit none -! ! arguments -! real(wp_), dimension(3), intent(in) :: xv0,anv0 -! real(wp_), dimension(3), intent(out) :: xvend -! real(wp_), intent(out) :: dstvac -! integer, intent(out) :: ivac -! ! local variables -! integer :: i -! real(wp_) :: st,rrm,zzm,smax -! real(wp_), dimension(3) :: walln -! logical :: plfound -! -! ! ivac=1 plasma hit before wall reflection -! ! ivac=2 wall hit before plasma -! ! ivac=-1 vessel (and thus plasma) never crossed -! -! call inters_linewall(xv0/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & -! nlim,smax,walln) -! smax=smax*1.0e2_wp_ -! rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) -! zzm=1.0e-2_wp_*xv0(3) -! if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then -! ! first wall interface is outside-inside -! if (dot_product(walln,walln) set_globals implicit none - type(antctrl_type) :: antp - type(eqparam_type) :: eqp - type(prfparam_type) :: prfp - type(outparam_type) :: outp - type(rtrparam_type) :: rtrp - type(hcdparam_type) :: hcdp - real(wp_), dimension(:), allocatable :: psrad, terad, derad, zfc - real(wp_), dimension(:), allocatable :: rv, zv, psinr, fpol, qpsi - real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim - real(wp_), dimension(:,:), allocatable :: psin - real(wp_) :: psia, rvac, rax, zax - real(wp_) :: fghz - real(wp_) :: x0, y0, z0, w1, w2, ri1, ri2, phiw, phir + ! gray_main subroutine arguments + type(gray_parameters) :: params ! Inputs + type(gray_data) :: data ! + type(gray_results) :: results ! Outputs + integer :: error ! Exit code - real(wp_) :: pec,icd + logical :: sum_mode = .false. - integer :: ierr - real(wp_), dimension(:), allocatable :: xrad, rhot, dpdv, jcd - real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx + ! Load the parameters and also copy them into + ! global variables exported by the gray_params + call read_parameters('gray_params.data', params) + call params_set_globals(params) - logical :: sum_mode = .FALSE. + ! Read the input data into set the global variables + ! of the respective module. Note: order matters. + call init_equilibrium(params, data) + call init_profiles(params, data) + call init_antenna(params%antenna) + call init_misc(params, data) - ! ------- sum mode variables ------- - real(wp_), dimension(:), allocatable :: jphi, currins, pins, rtin, rpin - integer :: i,j,k,n,ngam,irt - character(len=255) :: fn,dumstr - real(wp_), dimension(5) :: f48v - real(wp_) :: gam,alp,bet, jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - ! ---------------------------------- - - call read_params('gray_params.data', rtrp, hcdp, antp, eqp, rwallm, prfp, outp) - - ! ======= read input data BEGIN ======= - - !------------ equilibrium ------------ - if(eqp%iequil<2) then - call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim) - ! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) - psia = sign(one,qpsi(2)*fpol(1)) - else - call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & - rax,zax, rbnd,zbnd, rlim,zlim, & - eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) - call change_cocos(psia, fpol, qpsi, eqp%icocos, 3) - end if - ! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output - call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb) - ! ??? analytical only? change for numerical! - ! qpsi(1) = sign(qpsi(1),qpsi(1)*qpsi(2)*psia*fpol(1)) - ! qpsi(2) = sign(qpsi(2),psia*fpol(1)) - ! ---------------------------------- - - !------------- profiles ------------- - if(prfp%iprof == 0) then - call read_profiles_an(prfp%filenm, terad, derad, zfc) - else - call read_profiles(prfp%filenm, xrad, terad, derad, zfc) - allocate(psrad(size(xrad))) - if(prfp%irho == 0) then ! xrad==rhot - allocate(rhot(size(psinr))) - call setqphi_num(psinr,qpsi,psia,rhot) - call set_rhospl(sqrt(psinr),rhot) - deallocate(rhot) - psrad = frhopolv(xrad)**2 - else if(prfp%irho == 1) then ! xrad==rhop - psrad = xrad**2 - else - psrad = xrad - end if - deallocate(xrad) - end if - ! re-scale input data - call tene_scal(terad, derad, prfp%factte, prfp%factne, & - eqp%factb, prfp%iscal, prfp%iprof) - ! ---------------------------------- - - !------------- antenna -------------- - ! interpolate beam table if antctrl%ibeam>0 - select case (antp%ibeam) - case (2) - ! to be completed: now 1st beamd always selected, iox read from table - call read_beam2(antp%filenm, 1, antp%alpha, antp%beta, & - fghz, antp%iox, x0, y0, z0, & - w1, w2, ri1, ri2, phiw, phir) - case (1) - call read_beam1(antp%filenm, antp%alpha, antp%beta, & - fghz, x0, y0, z0, & - w1, w2, ri1, ri2, phiw, phir) - case default - call read_beam0(antp%filenm, & - fghz, x0, y0, z0, & - w1, w2, ri1, ri2, phiw, phir) - end select - ! ---------------------------------- - - !--------------- wall --------------- - ! set simple limiter if not read from EQDSK - ! need to clean up... - r0m=sqrt(x0**2+y0**2)*0.01_wp_ - dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_ - z0m=z0*0.01_wp_ - if (.not.allocated(rlim).or.rtrp%ipass<0) then - if (allocated(rlim)) deallocate(rlim) - if (allocated(zlim)) deallocate(zlim) - allocate(rlim(5)) - allocate(zlim(5)) - if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass - if(eqp%iequil<2) then - rmxm=(rv(1)+rv(2))*0.01_wp_ - else - rmxm=rv(size(rv)) - end if - call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim) - end if - ! ---------------------------------- - - ! ======= read input data END ======= - - ! ========================= MAIN SUBROUTINE CALL ========================= - allocate(dpdv(outp%nrho), jcd(outp%nrho)) if (sum_mode) then - allocate(jphi(outp%nrho), currins(outp%nrho), & - pins(outp%nrho), rtin(outp%nrho), rpin(outp%nrho)) + sum: block + real(wp_) :: pabs, icd, pec + real(wp_), dimension(:), allocatable :: dpdv, jcd, jphi + real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin + integer :: i, j, k, n, ngam, irt + character(len=255) :: filename + real(wp_), dimension(5) :: f48v + real(wp_) :: gam,alp,bet, jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx + allocate(jphi(params%output%nrho), currins(params%output%nrho), & + pins(params%output%nrho), rtin(params%output%nrho), & + rpin(params%output%nrho)) - open(100,file='filelist.txt',action='read',status='old') - read(100,*) n,ngam - do i=1,n - read(100,*) fn - open(100+i,file=fn,action='read',status='old') - do j=1,22 - read(100+i,*) dumstr - end do - end do - close(100) - - open(100+n+1,file='f48sum.txt',action='write',status='unknown') - open(100+n+2,file='f7sum.txt',action='write',status='unknown') - - do k=1,ngam - jphi=0.0_wp_ - jcd=0.0_wp_ - dpdv=0.0_wp_ - currins=0.0_wp_ - pins=0.0_wp_ - do j=1,outp%nrho - do i=1,n - read(100+i,*) gam,alp,bet,rpin(j),rtin(j),f48v(1:5),irt - jphi(j)=jphi(j)+f48v(1) - jcd(j)=jcd(j)+f48v(2) - dpdv(j)=dpdv(j)+f48v(3) - currins(j)=currins(j)+f48v(4) - pins(j)=pins(j)+f48v(5) + open(100, file='filelist.txt', action='read', status='old') + read(100, *) n, ngam + do i=1,n + read(100, *) filename + open(100 + i, file=filename, action='read', status='old') + do j=1,22 + read(100 + i, *) end do - write(100+n+1,'(10(1x,e16.8e3),i5)') gam,alp,bet,rpin(j),rtin(j), & - jphi(j),jcd(j),dpdv(j),currins(j),pins(j),irt end do - pec=pins(outp%nrho) - icd=currins(outp%nrho) - write(100+n+1,*) - call sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & - psrad,terad,derad,zfc,prfp, rlim,zlim, & - antp%beta,fghz,antp%alpha,antp%beta, & - (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, & - antp%psi, antp%chi, jphi,jcd,dpdv,currins,pins,pec,icd, & - jphip,dpdvp, rhotj,rhotjava,rhotp,rhotpav, & - drhotjava,drhotpav, ratjamx,ratjbmx, outp,rtrp,hcdp,ierr) - write(100+n+2,'(15(1x,e12.5),i5,4(1x,e12.5))') gam,alp,bet,icd,pec, & - jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - end do - do i=1,n+2 - close(100+i) - end do + close(100) + + open(100 + n+1, file='f48sum.txt', action='write', status='unknown') + open(100 + n+2, file='f7sum.txt', action='write', status='unknown') + + do k=1,ngam + jphi = zero + jcd = zero + dpdv = zero + currins = zero + pins = zero + do j=1,params%output%nrho + do i=1,n + read(100+i, *) gam, alp, bet, rpin(j), rtin(j), f48v(1:5), irt + jphi(j) = f48v(1) + jphi(j) + jcd(j) = f48v(2) + jcd(j) + dpdv(j) = f48v(3) + dpdv(j) + currins(j) = f48v(4) + currins(j) + pins(j) = f48v(5) + pins(j) + end do + write(100 + n+1,'(10(1x,e16.8e3),i5)') & + gam, alp, bet, rpin(j), rtin(j), & + jphi(j), jcd(j), dpdv(j), currins(j), pins(j), irt + end do + pec = pins(params%output%nrho) + icd = currins(params%output%nrho) + write(100 + n+1, *) + call sum_profiles(params, jphi, jcd, dpdv, currins, & + pins, pabs, icd, jphip, dpdvp, rhotj, & + rhotjava, rhotp, rhotpav, drhotjava, & + drhotpav, ratjamx, ratjbmx) + write(100 + n+2, '(15(1x,e12.5),i5,4(1x,e12.5))') & + gam, alp, bet, icd, pabs, jphip, dpdvp, & + rhotj, rhotjava, rhotp, rhotpav, & + drhotjava, drhotpav, ratjamx, ratjbmx + end do + do i=1,n+2 + close(100 + i) + end do + deallocate(dpdv, jcd, jphi, currins, pins, rtin, rpin) + end block sum else - call gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & - psrad,terad,derad,zfc,prfp, rlim,zlim, & - antp%power,fghz,antp%alpha,antp%beta, & - (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, & - antp%psi,antp%chi, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr) + call gray_main(params, data, results, error) end if - ! ======================================================================== - ! ======= control prints BEGIN ======= - if(ierr/=0) print*,' IERR = ', ierr - print*,' ' - print*,'Pabs (MW) = ', pec - print*,'Icd (kA) = ', icd*1.0e3_wp_ - ! ======= control prints END ======= + print '(a)' + print '(a,f9.4)', 'Pabs (MW)=', results%pabs + print '(a,f9.4)', 'Icd (kA)=', results%icd * 1.0e3_wp_ - ! ======= free memory BEGIN ======= - if(allocated(psrad)) deallocate(psrad) - if(allocated(terad)) deallocate(terad, derad, zfc) - if(allocated(rv)) deallocate(rv, zv, fpol, qpsi) - if(allocated(psin)) deallocate(psin, psinr) - if(allocated(rbnd)) deallocate(rbnd, zbnd) - if(allocated(rlim)) deallocate(rlim, zlim) - if(allocated(dpdv)) deallocate(dpdv, jcd) - if(sum_mode) deallocate(jphi, currins, pins, rtin, rpin) - ! ======= free memory END ====== -end program main_std + ! Free memory + call deinit_equilibrium(data%equilibrium) + call deinit_profiles(data%profiles) + call deinit_misc + deallocate(results%dpdv, results%jcd) + +contains + + subroutine init_equilibrium(params, data) + ! Reads the MHD equilibrium file (either in the G-EQDSK format + ! 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 + + implicit none + + ! subroutine arguments + type(gray_parameters), intent(inout) :: params + type(gray_data), intent(out) :: data + + 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) + + ! Set psia sign to give the correct sign to Iphi + ! (COCOS=3: psia<0 for Iphi>0) + data%equilibrium%psia = sign(one, data%equilibrium%qpsi(2) & + * data%equilibrium%fpol(1)) + else + ! Numerical equilibrium + call read_eqdsk(params%equilibrium, data%equilibrium) + call change_cocos(data%equilibrium, params%equilibrium%icocos, 3) + end if + + ! Rescale B, I and/or force their signs + 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)) + else + call set_eqspl(params%equilibrium, data%equilibrium) + end if + end subroutine init_equilibrium + + + 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 + + implicit none + + ! subroutine arguments + type(equilibrium_data), intent(inout) :: data + + ! Free the MHD equilibrium arrays + if (allocated(data%rv)) deallocate(data%rv, data%zv, data%fpol, data%qpsi) + if (allocated(data%psin)) deallocate(data%psin, data%psinr) + if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd) + if (allocated(data%rlim)) deallocate(data%rlim, data%zlim) + + ! Unset global variables of the `equilibrium` module + call unset_eqspl + call unset_rhospl + call unset_q + end subroutine deinit_equilibrium + + + subroutine init_profiles(params, data) + ! Reads the plasma kinetic profiles file (containing the elecron + ! temperature, density and plasma effective charge) and initialises + ! the respective GRAY data structure. + use gray_params, only : profiles_parameters, profiles_data + use equilibrium, only : frhopolv + use coreprofiles, only : read_profiles_an, read_profiles, tene_scal, & + set_prfan, set_prfspl + implicit none + + ! subroutine arguments + type(gray_parameters), intent(in) :: params + type(gray_data), intent(inout), target :: data + + ! local variables + type(profiles_parameters) :: profp + type(profiles_data), pointer :: profd + + ! Radial coordinate (depending on profp%irho: ρ_t, ρ_p, or ψ) + real(wp_), allocatable :: xrad(:) + + profp = params%profiles + profd => data%profiles + + if(params%profiles%iprof == 0) then + ! Analytical profiles + ! TODO: rewrite using derived type + call read_profiles_an(profp%filenm, profd%terad, profd%derad, profd%zfc) + else + ! Numerical profiles + call read_profiles(profp%filenm, xrad, profd%terad, profd%derad, profd%zfc) + + allocate(profd%psrad(size(xrad))) + + select case (profp%irho) + case (0) ! xrad is rhot + profd%psrad = frhopolv(xrad)**2 + case (1) ! xrad is rhop + profd%psrad = xrad**2 + case default ! xrad is psi + profd%psrad = xrad + end select + deallocate(xrad) + + end if + ! Rescale input data + ! TODO: rewrite using derived type + call tene_scal(profd%terad, profd%derad, profp%factte, profp%factne, & + params%equilibrium%factb, profp%iscal, profp%iprof) + + ! Set global variables + ! TODO: rewrite using derived type + if(params%profiles%iprof == 0) then + call set_prfan(profd%terad, profd%derad, profd%zfc) + else + call set_prfspl(profd%psrad, profd%terad, profd%derad, profd%zfc, & + profp%sspld, profp%psnbnd) + end if + end subroutine init_profiles + + + subroutine deinit_profiles(data) + ! Free all memory allocated by the init_profiles subroutine. + use gray_params, only : profiles_data + use coreprofiles, only : unset_prfspl + + implicit none + + ! subroutine arguments + type(profiles_data), intent(inout) :: data + + ! Free the plasma kinetic profiles arrays + if (allocated(data%psrad)) deallocate(data%psrad) + if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc) + + ! Unset global variables of the `coreprofiles` module + call unset_prfspl + end subroutine deinit_profiles + + + subroutine init_antenna(params) + ! Reads the wave launcher file (containing the wave frequency, launcher + ! position, direction and beam description) and initialises the respective + ! GRAY parameters. + use beams, only : read_beam0, read_beam1, read_beam2 + use gray_params, only : antenna_parameters + + implicit none + + ! subroutine arguments + type(antenna_parameters), intent(inout) :: params + + ! Note: α, β are loaded from gray_params.data + select case (params%ibeam) + case (2) + ! 2 degrees of freedom + ! w(z, α, β), 1/R(z, α, β) + ! FIXME: 1st beam is always selected, iox read from table + call read_beam2(params, beamid=1) + case (1) + ! 1 degree of freedom + ! w(z, α), 1/R(z, α) + call read_beam1(params) + case default + ! fixed w(z), 1/R(z) + call read_beam0(params) + end select + end subroutine init_antenna + + + subroutine init_misc(params, data) + ! Performs miscellanous initial tasks, before the gray_main subroutine. + use reflections, only : range2rect + use limiter, only : limiter_set_globals=>set_globals + + implicit none + + ! subroutine arguments + type(gray_parameters), intent(inout) :: params + type(gray_data), intent(inout) :: data + + ! Build a basic limiter when one is not provided by the EQDSK + if (.not. allocated(data%equilibrium%rlim) & + .or. params%raytracing%ipass < 0) then + block + real(wp_) :: rmxm, r0m, z0m, dzmx + r0m = sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2)* 0.01_wp_ + dzmx = abs(params%raytracing%ipass) * & + params%raytracing%dst * params%raytracing%nstep * 0.01_wp_ + z0m = params%antenna%pos(3) * 0.01_wp_ + + allocate(data%equilibrium%rlim(5)) + allocate(data%equilibrium%zlim(5)) + params%raytracing%ipass = abs(params%raytracing%ipass) + 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)) + end if + call range2rect(params%misc%rwall, max(r0m, rmxm), & + z0m - dzmx, z0m + dzmx, & + data%equilibrium%rlim, data%equilibrium%zlim) + end block + end if + + ! Set the global variables of the `limiter` module + call limiter_set_globals(data%equilibrium) + end subroutine init_misc + + + subroutine deinit_misc + ! Free all memory allocated by the init_misc subroutine. + use limiter, only : limiter_unset_globals=>unset_globals + + implicit none + + ! Unset the global variables of the `limiter` module + call limiter_unset_globals + end subroutine deinit_misc + + + subroutine sum_profiles(params, jphi, jcd, dpdv, currins, pins, pabs, icd, & + 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 dispersion, only : expinit + use gray_params, only : gray_parameters, print_parameters, & + headw, headl + use beams, only : launchangles2n, xgygcoeff + use magsurf_data, only : flux_average, dealloc_surfvec + use beamdata, only : init_btr, dealloc_beam + use pec, only : pec_init, postproc_profiles, dealloc_pec, & + rhop_tab, rhot_tab + use graycore, only : print_headers, print_finals, print_pec, & + print_bres, print_prof, print_maps, & + print_surfq + implicit none + + ! subroutine arguments + type(gray_parameters), intent(in) :: params + real(wp_), intent(in) :: pabs, icd + real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins + real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava, & + rhotp, rhotpav, drhotjava, drhotpav, & + ratjamx,ratjbmx + + ! local variables + real(wp_) :: ak0, bres, xgcn + real(wp_) :: chipol, psipol, st + real(wp_) :: drhotp, drhotj, dpdvmx, jphimx + + real(wp_), dimension(3) :: anv0 + real(wp_), dimension(:, :), pointer :: yw=>null(), ypw=>null(), gri=>null() + real(wp_), dimension(:, :, :), pointer :: xc=>null(), du1=>null(), ggri=>null() + + real(wp_), dimension(:, :), pointer :: psjki=>null(), ppabs=>null(), ccci=>null() + real(wp_), dimension(:), pointer :: tau0=>null(), alphaabs0=>null(), & + dids0=>null(), ccci0=>null() + real(wp_), dimension(:), pointer :: p0jk=>null() + complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() + integer, dimension(:), pointer :: iiv=>null() + + ! ======== set environment BEGIN ======== + ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) + call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) + + ! Compute the initial cartesian wavevector (anv0) + call launchangles2n(params%antenna, anv0) + + ! Initialise the ray variables (beamtracing) + call init_btr(params%raytracing, yw, ypw, xc, du1, & + gri, ggri, psjki, ppabs, ccci, & + tau0, alphaabs0, dids0, ccci0, & + p0jk, ext, eyt, iiv) + + ! Initialise the dispersion module + if(params%ecrh_cd%iwarm > 1) call expinit + + ! Initialise the magsurf_data module + call flux_average ! requires frhotor for dadrhot,dvdrhot + + ! Initialise the output profiles + call pec_init(params%output%ipec) + ! ======= set environment END ====== + + ! ======== pre-proc prints BEGIN ======== + block + ! Parameters log in file headers + character(len=headw), dimension(headl) :: strheader + call print_parameters(params, strheader) + call print_headers(strheader) + end block + + ! Print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout + call print_surfq([1.5_wp_, 2.0_wp_]) + + ! Print ne, Te, q, Jphi versus psi, rhop, rhot + call print_bres(bres) + call print_prof + call print_maps(bres, xgcn, & + 0.01_wp_*sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2), & + sin(params%antenna%beta*degree)) + ! ========= pre-proc prints END ========= + + ! Print power and current density profiles + call print_pec(rhop_tab, rhot_tab, jphi, jcd, & + dpdv, currins, pins, index_rt=1) + ! Compute profiles width + call postproc_profiles(pabs, icd, rhot_tab, dpdv, jphi, & + rhotpav, drhotpav, rhotjava, drhotjava, & + dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, & + dpdvmx, jphimx, ratjamx, ratjbmx) + ! Print 0D results + call print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, & + drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, & + drhotj, ratjamx, ratjbmx, st, psipol, chipol, & + 1, params%antenna%power, cpl1=zero, cpl2=zero) + + ! Free memory + call dealloc_surfvec ! for fluxval + call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, & + tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv) + call dealloc_pec + end subroutine sum_profiles + +end program main