diff --git a/Makefile.single b/Makefile.single new file mode 100644 index 0000000..242ade8 --- /dev/null +++ b/Makefile.single @@ -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) $< + diff --git a/src/gray-externals.f b/src/gray-externals.f index 8197a4f..71369b6 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -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 diff --git a/src/gray-single.f90 b/src/gray-single.f90 new file mode 100644 index 0000000..eee787f --- /dev/null +++ b/src/gray-single.f90 @@ -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 diff --git a/src/gray.f b/src/gray.f index 58fe4ee..5c2390a 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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