Fixed memory allocation errors in read_beams. Added driver for single beam simulations, with beam parameters read from external file.
This commit is contained in:
parent
f997ec1eb3
commit
b9b687be88
165
Makefile.single
Normal file
165
Makefile.single
Normal file
@ -0,0 +1,165 @@
|
||||
# Executable name
|
||||
EXE=gray-single
|
||||
|
||||
# Objects list
|
||||
OBJ= gray-single.o gray_main.o gray-externals.o grayl.o beamdata.o \
|
||||
const_and_precisions.o green_func_p.o reflections.o scatterspl.o
|
||||
|
||||
# Alternative search paths
|
||||
vpath %.f90 src
|
||||
vpath %.f src
|
||||
|
||||
########################################
|
||||
# Set up compiler specific environment #
|
||||
########################################
|
||||
# set default compiler
|
||||
# --------------------
|
||||
FC=gfortran
|
||||
|
||||
# Set compiler debug and optimization flags
|
||||
# -----------------------------------------
|
||||
FFLAGS= -O0 -finit-real=nan -fcheck=all
|
||||
|
||||
#####################################
|
||||
# Set up RULES specific environment #
|
||||
#####################################
|
||||
|
||||
DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
|
||||
|
||||
# suffixes
|
||||
# --------
|
||||
.SUFFIXES: .f90 .f .o .mod
|
||||
|
||||
# PHONY targets
|
||||
# -------------
|
||||
.PHONY: clean install
|
||||
|
||||
all: $(EXE)
|
||||
|
||||
# Build executable from object files
|
||||
$(EXE): $(OBJ)
|
||||
$(FC) $(FFLAGS) -o $@ $^
|
||||
|
||||
# create rule for compiling f90 files
|
||||
# -----------------------------------
|
||||
%.o: %.f90
|
||||
$(FC) -c $(FFLAGS) $<
|
||||
|
||||
# create rule for compiling f files
|
||||
# ---------------------------------
|
||||
%.o: %.f
|
||||
$(FC) -c $(FFLAGS) $<
|
||||
|
||||
# create rule for files requiring preprocessing
|
||||
gray-externals.o: gray-externals.f
|
||||
$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $<
|
||||
|
||||
# make 'clean' option
|
||||
# -------------------
|
||||
clean:
|
||||
rm -rf *.o *.mod $(EXE)
|
||||
|
||||
# Dependencies
|
||||
# ------------
|
||||
gray-single.o: gray_main.o grayl.o
|
||||
gray_main.o: const_and_precisions.o
|
||||
gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
|
||||
green_func_p.o: const_and_precisions.o
|
||||
scatterspl.o: const_and_precisions.o
|
||||
beamdata.o: const_and_precisions.o
|
||||
|
||||
## library name
|
||||
## ------------
|
||||
#LIBNAME=$(JLIBDIR)/libgray.a
|
||||
#
|
||||
## List source and object files
|
||||
## ----------------------------
|
||||
#FSRC=$(wildcard *.f *.f90)
|
||||
#FOBJ0=$(FSRC:.f=.o)
|
||||
#FOBJ=$(FOBJ0:.f90=.o)
|
||||
#
|
||||
#########################################
|
||||
## Set up compiler specific environment #
|
||||
#########################################
|
||||
## set default compiler
|
||||
## --------------------
|
||||
#FC=$(F90)
|
||||
#
|
||||
## Set compiler debug and optimization flags
|
||||
## -----------------------------------------
|
||||
#ifeq ("$(DBG)","")
|
||||
# FFLAGS= $(AUTO) -O3
|
||||
#else
|
||||
# FFLAGS= $(AUTO) -O0 -g -Minform=inform -Mbounds -Mchkptr
|
||||
#endif
|
||||
#
|
||||
## Set include directories
|
||||
## -----------------------
|
||||
##INC=-I$(JCOMMON_STD)
|
||||
#INC=
|
||||
#
|
||||
######################################
|
||||
## Set up RULES specific environment #
|
||||
######################################
|
||||
#
|
||||
## suffixes
|
||||
## --------
|
||||
#.SUFFIXES: .f90 .f .o .mod
|
||||
#
|
||||
## PHONY targets
|
||||
## -------------
|
||||
#.PHONY: all clean realclean $(MASTER_RULES) gray
|
||||
#
|
||||
## default target
|
||||
## --------------
|
||||
#all:
|
||||
# @($(MAKE) -C ../ all) || exit $$?
|
||||
#
|
||||
## Set rules passed to top-level Makefile
|
||||
## --------------------------------------
|
||||
#$(MASTER_RULES):
|
||||
# @($(MAKE) -C ../ $@) || exit $$?
|
||||
#
|
||||
## rule for libgray.a
|
||||
## --------------------
|
||||
#$(LIBNAME): $(FOBJ)
|
||||
# ar vr $(LIBNAME) $?
|
||||
#
|
||||
## synonym for libgray.a
|
||||
## -----------------------
|
||||
#gray: $(LIBNAME)
|
||||
#
|
||||
## create rule for compiling f90 files
|
||||
## -----------------------------------
|
||||
#%.o: %.f90
|
||||
# $(FC) -c $(FFLAGS) $(INC) $<
|
||||
#
|
||||
## create rule for compiling f files
|
||||
## ---------------------------------
|
||||
#%.o: %.f
|
||||
# $(FC) -c $(FFLAGS) $(INC) $<
|
||||
#
|
||||
## make 'realclean' option
|
||||
## -----------------------
|
||||
#realclean:
|
||||
# rm -f *.o *.mod $(LIBNAME)
|
||||
#
|
||||
## make 'clean' option
|
||||
## -------------------
|
||||
#clean:
|
||||
# rm -f *.o
|
||||
#
|
||||
## Dependencies
|
||||
## ------------
|
||||
#gray_main.o: const_and_precisions.o
|
||||
#gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
|
||||
#green_func_p.o: const_and_precisions.o
|
||||
#scatterspl.o: const_and_precisions.o
|
||||
#beamdata.o: const_and_precisions.o
|
||||
#
|
||||
## Special rule to preprocess source file and include svn revision
|
||||
## ---------------------------------------------------------------
|
||||
#DIRECTIVES = -DREVISION="'$(shell svnversion)'"
|
||||
#gray-externals.o:gray-externals.f
|
||||
# $(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC) $<
|
||||
|
||||
@ -533,7 +533,7 @@ c angles read from values passed as arguments (in rad)
|
||||
alpha0 = alphain
|
||||
beta0 = betain
|
||||
read(602,*) fghz
|
||||
read(602,*) p0mw
|
||||
read(602,*) dummy !p0mw
|
||||
c power value overwritten with value passed as argument (in W)
|
||||
p0mw=powin*1.d-6
|
||||
c
|
||||
@ -560,6 +560,7 @@ c ibeam=1 :read data from file simple astigmatic beam
|
||||
c ibeam=2 :read data from file general astigmatic beam
|
||||
read(602,*) ibeam
|
||||
read(602,*) filenmbm
|
||||
filenmbm='graybeam.data'
|
||||
c
|
||||
c iequil=0 :vacuum
|
||||
c iequil=1 :analytical equilibrium
|
||||
@ -691,7 +692,6 @@ c
|
||||
c read data for beam from file if ibeam>0
|
||||
c
|
||||
if(ibeam.gt.0) then
|
||||
filenmbm='graybeam.data'
|
||||
call read_beams(beamid,iox)
|
||||
else
|
||||
zrcsi=0.5d0*ak0*w0csi**2
|
||||
@ -937,8 +937,7 @@ c ibeam=1 simple astigmatic beam
|
||||
c ibeam=2 general astigmatic beam
|
||||
c
|
||||
nfbeam=603
|
||||
lbm=length(filenmbm)
|
||||
open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam)
|
||||
open(file=trim(filenmbm),status= 'old',unit=nfbeam)
|
||||
c
|
||||
c#######################################################################################
|
||||
if(ibeam.eq.1) then
|
||||
@ -947,7 +946,8 @@ c###############################################################################
|
||||
nbeta = 1
|
||||
allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer),
|
||||
. waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer),
|
||||
. phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer))
|
||||
. phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer),
|
||||
. xcoord(nisteer),ycoord(nisteer))
|
||||
c
|
||||
c c====================================================================================
|
||||
do i=1,nisteer
|
||||
@ -1014,7 +1014,8 @@ c # of elements in beam data grid
|
||||
c
|
||||
allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer),
|
||||
. waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer),
|
||||
. phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer))
|
||||
. phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer),
|
||||
. xcoord(nisteer),ycoord(nisteer))
|
||||
c
|
||||
c c====================================================================================
|
||||
c beam data grid reading
|
||||
@ -1128,7 +1129,7 @@ c
|
||||
. txrci1,crci1,fp,wrk,lwrk,iwrk,ier)
|
||||
c
|
||||
call dierckx_curfit(iopt,nxcoord,xcoord,rci2v,w,
|
||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci21,
|
||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2,
|
||||
. txrci2,crci2,fp,wrk,lwrk,iwrk,ier)
|
||||
c
|
||||
call dierckx_curfit(iopt,nxcoord,xcoord,phi1v,w,
|
||||
@ -1255,8 +1256,8 @@ c check if (xcoord0, ycoord0) is out of table range
|
||||
c incheck = 1 inside / 0 outside
|
||||
if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0))
|
||||
. incheck = 1
|
||||
deallocate(wrk,iwrk)
|
||||
end if
|
||||
deallocate(wrk,iwrk)
|
||||
c#######################################################################################
|
||||
c
|
||||
c
|
||||
@ -1267,7 +1268,7 @@ c c============================================================================
|
||||
call dierckx_splev(txycoord,nxycoord,cycoord,kx,xcoord0,
|
||||
. ycoord0,1,ier)
|
||||
c
|
||||
call dierckx_splev(txycoord,nxycoord,cwaist1,kx,xcoord0,wcsi,
|
||||
call dierckx_splev(txwaist1,nxwaist1,cwaist1,kx,xcoord0,wcsi,
|
||||
. 1,ier)
|
||||
c
|
||||
call dierckx_splev(txwaist2,nxwaist2,cwaist2,kx,xcoord0,weta,
|
||||
@ -1446,6 +1447,8 @@ c print new (alpha0, beta0)
|
||||
write(*,10) ' new (alpha0,beta0) = (', xcoord0,
|
||||
. ',', ycoord0, ')'
|
||||
10 format(a27,f10.5,a1,f10.5,a1)
|
||||
deallocate(xpolygA, ypolygA, xpolygC, ypolygC,
|
||||
. xpolygB, ypolygB, xpolygD, ypolygD)
|
||||
end if
|
||||
c c====================================================================================
|
||||
c
|
||||
@ -1483,6 +1486,7 @@ c
|
||||
c
|
||||
call dierckx_bispev(txz0,nxz0,tyz0,nyz0,cz0,kx,ky,
|
||||
. (/xcoord0/),1,(/ycoord0/),1,z00,wrk,lwrk,iwrk,kwrk,ier)
|
||||
deallocate(wrk,iwrk)
|
||||
c c----------------------------------------------------------------------------------
|
||||
else
|
||||
c c----------------------------------------------------------------------------------
|
||||
@ -1500,6 +1504,18 @@ c c============================================================================
|
||||
end if
|
||||
c#######################################################################################
|
||||
c
|
||||
if(fdeg.ne.0) then
|
||||
deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2,
|
||||
. txwaist2, crci1, txrci1, crci2, txrci2, cphi1, txphi1,
|
||||
. cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w)
|
||||
else
|
||||
deallocate(cwaist1, txwaist1, tywaist1,
|
||||
. cwaist2, txwaist2, tywaist2,
|
||||
. crci1, txrci1, tyrci1, crci2, txrci2, tyrci2,
|
||||
. cphi1, txphi1, typhi1, cphi2, txphi2, typhi2,
|
||||
. cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0,
|
||||
. wrk2, xpolyg, ypolyg, w)
|
||||
end if
|
||||
c
|
||||
c#######################################################################################
|
||||
c set correct values for alpha, beta
|
||||
@ -1511,6 +1527,8 @@ c set correct values for alpha, beta
|
||||
beta0 = ycoord0
|
||||
end if
|
||||
c#######################################################################################
|
||||
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v,
|
||||
. phi2v,x00v,y00v,z00v,xcoord,ycoord)
|
||||
c
|
||||
return
|
||||
end
|
||||
|
||||
306
src/gray-single.f90
Normal file
306
src/gray-single.f90
Normal file
@ -0,0 +1,306 @@
|
||||
|
||||
! (ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd,
|
||||
! . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne,
|
||||
! . zeff, qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec,
|
||||
! . icd, ierr)
|
||||
|
||||
implicit real*8 (a-h,o-z)
|
||||
parameter(pi=3.14159265358979d0)
|
||||
! input arguments
|
||||
integer :: ijetto, mr, mz, nbnd, nrho, nbeam
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: r, z, psinr
|
||||
real*8, DIMENSION(:,:), ALLOCATABLE :: psin, psi
|
||||
real*8 psiax, psibnd, rax, zax
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: rbnd, zbnd
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: psijet, f, qsf, te, dne, zeff
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: teitpl, dneitpl, zeffitpl, cte, ctdn, ctz
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: powin, alphain, betain
|
||||
! do loop arguments
|
||||
real*8 :: powloop,alphaloop,betaloop
|
||||
! do loop output arguments
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: dpdv, jcd
|
||||
real*8 :: pec, icd
|
||||
! gray_main output arguments
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: dpdvloop, jcdloop
|
||||
real*8 :: pecloop, icdloop
|
||||
integer :: ierr
|
||||
! unused read variales
|
||||
real*8, DIMENSION(:), ALLOCATABLE :: pres,ffprim,pprim,rlim,zlim
|
||||
character(LEN=48) :: stringa
|
||||
!!!!!!!!!! VARIABILI PER TEST !!!!!!!!!!
|
||||
real*8,dimension(2) :: alphalim, betalim
|
||||
real*8 :: delta
|
||||
integer :: alphanum, betanum, len
|
||||
character(LEN=48) :: profile, eqfile
|
||||
!!!!!!!!!! VARIABILI PER TEST !!!!!!!!!!
|
||||
!
|
||||
! local variables
|
||||
! real*8 fgray(nrho), qgray(nrho), jcdgry(nrho), icdgry
|
||||
!
|
||||
! input arguments
|
||||
!
|
||||
! ijetto Equilibrium source (1 EFIT, 2 ESCO)
|
||||
! If IJETTO=2, then PSIN values are valid only inside
|
||||
! plasma boudary (PSIN=0 outside)
|
||||
! mr Size of flux map grid in R direction
|
||||
! mz Size of flux map grid in Z direction
|
||||
! mrd Leading dimension of the psin(:,:) array, mrd>mr
|
||||
! r R coordinates of flux map grid points [m]
|
||||
! z Z coordinates of flux map grid points [m]
|
||||
! psin Normalised poloidal flux psin=(psi-psiax)/(psibnd-psiax)
|
||||
! on the (R, Z) grid.
|
||||
! psiax Poloidal flux on axis [Wb rad-1]
|
||||
! psibnd Poloidal flux on boundary [Wb rad-1]
|
||||
! rax R coordinate of magnetic axis [m]
|
||||
! zax Z coordinate of magnetic axis [m]
|
||||
! nbnd Number of points in plasma boundary contour
|
||||
! rbnd R coordinates of plasma boundary contour [m]
|
||||
! zbnd Z coordinates of plasma boundary contour [m]
|
||||
!
|
||||
! nrho Number of points in JETTO rho grid -
|
||||
! psijet Normalised poloidal flux on JETTO radial grid
|
||||
! f Poloidal current stream function f=Bphi*R on JETTO
|
||||
! radial grid [T m]
|
||||
! te Electron temperature on JETTO radial grid [eV]
|
||||
! dne Electron density on JETTO radial grid [m-3]
|
||||
! zeff Effective nuclear charge Zeff on JETTO radial grid
|
||||
! qsf Safety factor on JETTO radial grid
|
||||
!
|
||||
! nbeam number of beams injected
|
||||
! powin Input ECRH power array [W] (powin(i) =< 0 means i-th beam is unused)
|
||||
! alphin Beams poloidal injection angles array [rad]
|
||||
! betain Beams toroidal injection angles array [rad]
|
||||
!
|
||||
! output arguments
|
||||
!
|
||||
! dpdv Absorbed EC power density on JETTO radial grid [W m-3]
|
||||
! jcd EC driven flux averaged current density on JETTO
|
||||
! radial grid [A m-2]
|
||||
! pec Total absorbed EC power [W]
|
||||
! icd Total EC driven current [A]
|
||||
! ierr Return code. IERR>0 on error
|
||||
! ierr = 90-93: error computing integrals for current drive
|
||||
! ierr = 94: absorption coefficient alpha < 0
|
||||
! ierr = 97: parallel comp. refract. idx N//>0.99 (warning)
|
||||
! ierr = 98: parallel comp. refract. idx N//>1.05
|
||||
!
|
||||
! JETTO coordinate system assumes toroidal angle increasing CW
|
||||
! in GRAY toroidal angle increases CCW --> adapt signs on input data
|
||||
|
||||
!!!!!!!!!! VARIABILI PER TEST !!!!!!!!!!
|
||||
|
||||
allocate(powin(1),alphain(1),betain(1))
|
||||
nbeam=1
|
||||
|
||||
ijetto=1
|
||||
|
||||
! lettura gray.data
|
||||
open(99,file='gray.data')
|
||||
read(99,*) alphain(1),betain(1)
|
||||
read(99,*)
|
||||
read(99,*) powin(1)
|
||||
powin(1)=powin(1)*1.d6
|
||||
do i=1,13
|
||||
read(99,*)
|
||||
end do
|
||||
read(99,*) eqfile
|
||||
read(99,*) ipsinorm
|
||||
read(99,*) factb
|
||||
read(99,*) profile
|
||||
read(99,*)
|
||||
read(99,*) sgnbphi,sgniphi,icocos
|
||||
close(99)
|
||||
|
||||
! lettura profili di flusso poloidale, temperatura, densità e carica nucleare
|
||||
len=length(profile)
|
||||
open(98,file=profile(1:len)//'.prf',action='read')
|
||||
read(98,*) nrho
|
||||
allocate(psijet(nrho),te(nrho),dne(nrho),zeff(nrho))
|
||||
read(98,*) psijet0,te0,dne0,zeff0
|
||||
if(psijet0.ne.0.0d0) psijet0=0.0d0
|
||||
psijet(1)=psijet0
|
||||
te(1)=te0
|
||||
dne(1)=dne0
|
||||
zeff(1)=zeff0
|
||||
do i=2,nrho
|
||||
read(98,*) psijeti,tei,dnei,zeffi !,rhoi,qsfi
|
||||
psijet(i)=psijeti
|
||||
te(i)=tei
|
||||
dne(i)=dnei
|
||||
zeff(i)=zeffi
|
||||
end do
|
||||
close(98)
|
||||
|
||||
! lettura dati campo magnetico
|
||||
len=length(eqfile)
|
||||
open(99,file=eqfile(1:len)//'.eqdsk',action='read')
|
||||
|
||||
if(icocos.gt.0) then
|
||||
read (99,'(a48,3i4)') stringa,idum,mr,mz
|
||||
else
|
||||
read (99,*) mr,mz
|
||||
end if
|
||||
|
||||
allocate(r(mr),z(mz),psinr(mr),psi(mr,mz),psin(mr,mz),f(mr),qsf(mr))
|
||||
allocate(pres(mr),ffprim(mr),pprim(mr))
|
||||
allocate(dpdv(mr),jcd(mr),dpdvloop(mr), jcdloop(mr))
|
||||
|
||||
if(icocos.gt.0) then
|
||||
read (99,'(5e16.9)') drnr1,dznz1,rcen,rleft,zmid
|
||||
read (99,'(5e16.9)') rax,zax,psiax,psibnd,btrcen
|
||||
read (99,'(5e16.9)') current,xdum,xdum,xdum,xdum
|
||||
read (99,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
|
||||
read (99,'(5e16.9)') (f(i),i=1,mr)
|
||||
read (99,'(5e16.9)') (pres(i),i=1,mr)
|
||||
read (99,'(5e16.9)') (ffprim(i),i=1,mr)
|
||||
read (99,'(5e16.9)') (pprim(i),i=1,mr)
|
||||
if (ipsinorm.eq.0) then
|
||||
read (99,'(5e16.9)') ((psi(i,j),i=1,mr),j=1,mz)
|
||||
else
|
||||
read (99,'(5e16.9)') ((psin(i,j),i=1,mr),j=1,mz)
|
||||
end if
|
||||
read (99,'(5e16.9)') (qsf(i),i=1,mr)
|
||||
read (99,*) nbnd,nlim
|
||||
allocate(rbnd(nbnd),zbnd(nbnd),rlim(nlim),zlim(nlim))
|
||||
read (99,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd)
|
||||
read (99,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
|
||||
else
|
||||
read (99,*) drnr1,dznz1,rcen,rleft,zmid
|
||||
read (99,*) rax,zax,psiax,psibnd,btrcen
|
||||
read (99,*) current,xdum,xdum,xdum,xdum
|
||||
read (99,*) xdum,xdum,xdum,xdum,xdum
|
||||
read (99,*) (f(i),i=1,mr)
|
||||
read (99,*) (pres(i),i=1,mr)
|
||||
read (99,*) (ffprim(i),i=1,mr)
|
||||
read (99,*) (pprim(i),i=1,mr)
|
||||
if (ipsinorm.eq.0) then
|
||||
read (99,*) ((psi(i,j),i=1,mr),j=1,mz)
|
||||
else
|
||||
read (99,*) ((psin(i,j),i=1,mr),j=1,mz)
|
||||
end if
|
||||
read (99,*) (qsf(i),i=1,mr)
|
||||
read (99,*) nbnd,nlim
|
||||
allocate(rbnd(nbnd),zbnd(nbnd),rlim(nlim),zlim(nlim))
|
||||
read (99,*) (rbnd(i),zbnd(i),i=1,nbnd)
|
||||
read (99,*) (rlim(i),zlim(i),i=1,nlim)
|
||||
end if
|
||||
close(99)
|
||||
|
||||
! adeguamento alle convenzioni e normalizzazioni
|
||||
icocosmod=mod(icocos,10)
|
||||
|
||||
! icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
|
||||
if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
|
||||
btrcen=-btrcen
|
||||
current=-current
|
||||
do i=1,mr
|
||||
f(i)=-f(i)
|
||||
end do
|
||||
end if
|
||||
|
||||
if (icocosmod.ne.0) sgnbphi=sign(1.0d0,f(mr))
|
||||
btrcen=sign(btrcen,sgnbphi)
|
||||
|
||||
do i=1,mr
|
||||
f(i)=sign(f(i),sgnbphi)*factb
|
||||
end do
|
||||
|
||||
! icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
|
||||
! icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
|
||||
if (icocosmod.eq.1 .or. icocosmod.eq.4 .or.icocosmod.eq.5 .or. icocosmod.eq.8) then
|
||||
psibnd=-psibnd
|
||||
psiax=-psiax
|
||||
if (ipsinorm.eq.0) then
|
||||
do j=1,mz
|
||||
do i=1,mr
|
||||
psi(i,j)=-psi(i,j)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
|
||||
psia0=psibnd-psiax
|
||||
|
||||
! icocos=0: adapt psi to force Ip sign, otherwise maintain psi
|
||||
if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0)
|
||||
current=sign(current,sgniphi)
|
||||
psia=sign(psia0,-sgniphi)*factb
|
||||
|
||||
! icocos>10: input psi is in Wb
|
||||
! icocos<10: input psi is in Wb/rad (gray convention)
|
||||
if (icocos.ge.10) psia=psia/(2.0d0*pi)
|
||||
|
||||
if(ipsinorm.eq.0) then
|
||||
do j=1,mz
|
||||
do i=1,mr
|
||||
psin(i,j)=(psi(i,j)-psiax)/(psia0)
|
||||
psi(i,j)=psin(i,j)*psia
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,mz
|
||||
do i=1,mr
|
||||
psi(i,j)=psin(i,j)*psia
|
||||
enddo
|
||||
enddo
|
||||
end if
|
||||
|
||||
dr=drnr1/dble(mr-1)
|
||||
dz=dznz1/dble(mz-1)
|
||||
r(1)=rleft
|
||||
z(1)=zmid-dznz1/2.0d0
|
||||
dpsinr=1.0d0/dble(mr-1)
|
||||
|
||||
do i=1,mr
|
||||
psinr(i)=dpsinr*(i-1)
|
||||
r(i)=r(1)+(i-1)*dr
|
||||
end do
|
||||
|
||||
do j=1,mz
|
||||
z(j)=z(1)+(j-1)*dz
|
||||
end do
|
||||
|
||||
! interpolazione Te, dne, Zeff (psijet -> psinr, nrho -> mr elementi)
|
||||
|
||||
allocate(teitpl(mr),dneitpl(mr),zeffitpl(mr))
|
||||
allocate(cte(4*nrho),ctdn(4*nrho),ctz(4*nrho))
|
||||
|
||||
iopt=0
|
||||
call difcsg(psijet,te,nrho,iopt,cte,ier)
|
||||
call difcsg(psijet,dne,nrho,iopt,ctdn,ier)
|
||||
call difcsg(psijet,zeff,nrho,iopt,ctz,ier)
|
||||
|
||||
do i=1,mr
|
||||
call vlocate(psijet,nrho,psinr(i),k)
|
||||
k = max(1,min(k,nrho-1))
|
||||
dal = psinr(i)-psijet(k)
|
||||
|
||||
teitpl(i) = fspli(cte,nrho,k,dal)*1000.d0
|
||||
dneitpl(i) = fspli(ctdn,nrho,k,dal)*1.d19
|
||||
zeffitpl(i) = fspli(ctz,nrho,k,dal)
|
||||
end do
|
||||
|
||||
! inizializzazione variabili output
|
||||
dpdv = 0
|
||||
jcd = 0
|
||||
pec = 0
|
||||
icd = 0
|
||||
|
||||
do j=1,nbeam
|
||||
! potenza ECRH e angoli per fascio j-esimo
|
||||
powloop = powin(j)
|
||||
alphaloop = alphain(j)
|
||||
betaloop = betain(j)
|
||||
|
||||
call gray_main(ijetto, mr, mz, r, z, psin, 0.0d0, psia, &
|
||||
rax, zax, nbnd, rbnd, zbnd, mr, psinr, f, teitpl, dneitpl, &
|
||||
zeffitpl, qsf, j, powloop, alphaloop, betaloop, &
|
||||
dpdvloop, jcdloop, pecloop, icdloop, ierr)
|
||||
|
||||
dpdv = dpdv + dpdvloop
|
||||
jcd = jcd + jcdloop
|
||||
pec = pec + pecloop
|
||||
icd = icd + icdloop
|
||||
end do
|
||||
|
||||
end
|
||||
@ -106,6 +106,8 @@ c f(i)=-f(i)
|
||||
c qsf(i)=-qsf(i)
|
||||
end do
|
||||
icd=-icd
|
||||
c
|
||||
close(604,607,608,609,612,617,630,631,633,644,645,646,648)
|
||||
c
|
||||
return
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user