diff --git a/Makefile b/Makefile index db3080b..d525689 100644 --- a/Makefile +++ b/Makefile @@ -12,6 +12,8 @@ vpath %.f src FC=gfortran FFLAGS=-O3 +DIRECTIVES = -DREVISION="'$(shell svnversion)'" + all: $(EXE) # Build executable from object files @@ -28,7 +30,7 @@ itm_constants.o: itm_types.o $(FC) $(FFLAGS) -c $< gray.o:gray.f green_func_p.o - $(FC) $(FFLAGS) -c $< + $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< grayl.o:grayl.f $(FC) $(FFLAGS) -c $^ diff --git a/src/gray.f b/src/gray.f index d889291..7aaaf12 100644 --- a/src/gray.f +++ b/src/gray.f @@ -42,7 +42,7 @@ c second pass into plasma call prfile call vectinit2 call paraminit - call ic_rt0 + call ic_rt2 call gray_integration call after_gray_integration pabstott=pabstott+pabstot @@ -54,7 +54,7 @@ c second pass into plasma call prfile call vectinit2 call paraminit - call ic_rt0 + call ic_rt2 call gray_integration call after_gray_integration pabstott=pabstott+pabstot @@ -64,7 +64,7 @@ c second pass into plasma 999 continue print*,' ' print*,' IERR = ', ierr - print*,'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3 + write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3 stop end @@ -235,6 +235,7 @@ c ierr=0 psjki(j,k,i)=psinv rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 + if(j.eq.1) rrm11=rrm if (iwarm.gt.0.and.i.gt.1) then if(psinv.ge.0.and.psinv.le.1.0d0) then @@ -285,7 +286,7 @@ c print*,' ' else powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2) end if - print*,'Reflected power fraction =',powrfl + write(6,*) 'Reflected power fraction =',powrfl iov(j,k)=3 yyrfl(j,k,1)=xvrfl(1) yyrfl(j,k,2)=xvrfl(2) @@ -320,7 +321,7 @@ c print*,' ' psimin=psjki(1,1,i) if(nray.gt.1) psimin=min(psimin,minval(psjki(2:nray,1:ktx,i))) - if(psimin.gt.1.0d0.and.rrm.gt.rcen.and.index_rt.gt.1) + if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1) . istop=1 if(iovmin.eq.3) istop=1 @@ -492,15 +493,16 @@ c write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi'// .' N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' - write(8,*) ' #istep j k xt yt zt rt ' - write(9,*) ' #istep j k xt yt zt rt ' - write(17,*) ' #sst ddr_11 ddr_Nr1 ddi_Nr1' + write(8,*) ' #istep j k xt yt zt rt psin' + write(9,*) ' #istep j k xt yt zt rt psin' + write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha' - write(12,*) ' #i sst rhop w1 w2 w1a w2a' + write(12,*) ' #i sst psi w1 w2' + write(7,*)'#Icd Pa Jphimx dPdVmx'// + .'rhotj rhotjava rhotp rhotpav '// + .'drhotjava drhotpav stmx polpsi polchi index_rt' + write(48,*) '#psi rhot Jphi Jcda dPdV Icdins Pins P% index_rt' -c write(29,*) '#beta alfa qqom uuom vvom psiom chiom'// -c .' qqxm uuxm vvxm psixm chixm' -c write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm' else If(index_rt.eq.3) then @@ -509,6 +511,7 @@ c write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm' write(9,*) ' ' write(17,*) ' ' write(12,*) ' ' + write(48,*) ' ' end if end if @@ -521,6 +524,8 @@ c use green_func_p, only:Setup_SpitzFunc implicit real*8 (a-h,o-z) real*8 me + character*8 wdat + character*10 wtim character*24 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) @@ -816,47 +821,58 @@ c versus psi, rhop, rhot if (iequil.eq.1) call surf_anal - if (iequil.eq.2.and.iprof.eq.1) then nfil=78 - open(nfil,file='tit.txt',status= 'unknown') - write(nfil,905) filenmeqq - write(nfil,907) filenmprf - write(nfil,911) fghz - write(nfil,914) btrcen, sgnbphi,sgniphi,icocos - write(nfil,900) nray, ktx, rmx + + open(nfil,file='headers.txt',status= 'unknown') + call date_and_time(wdat,wtim) + write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8), + . wtim(1:2),wtim(3:4),wtim(5:6) + write(nfil,904) REVISION + if (iequil.eq.2) then + write(nfil,905) filenmeqq + else + if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq + if (iequil.eq.0) write(nfil,'("# VACUUM CASE")') + end if + if (iprof.eq.1) then + write(nfil,907) filenmprf + else + write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff + end if + write(nfil,911) fghz,p0mw,iox write(nfil,902) x00,y00,z00 + write(nfil,908) alfac,betac if(ibeam.ge.1) write(nfil,909) filenmbm if(ibeam.eq.0) write(nfil,903) w0xt,w0yt,pw0xt,pw0yt,awr + write(nfil,900) nray, ktx, rmx write(nfil,901) igrad,iwarm,ilarm,ieccd,idst - if(ieccd.eq.1) write(nfil,912) - if(ieccd.eq.11) write(nfil,913) - write(nfil,904) p0mw,iox + write(nfil,915) dst,nstep + write(nfil,914) sgnbphi,sgniphi,icocos write(nfil,906) factb,factt,factn,iscal - write(nfil,910) ipec,nnd,ipsinorm,sspl,psdbnd -c write(nfil,915) psi15,sqrt(psi15),rhot15 -c write(nfil,920) psi2,sqrt(psi2),rhot2 + write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm write(nfil,*) ' ' close(nfil) - end if return -900 format('# Nray ktx rmx',2i5,1x,e12.5) -901 format('# igrad iwarm ilarm ieccd idst',6i5) -902 format('# X Y Z launching point (cm) : ',3(1x,e12.5)) -903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,e12.5)) -904 format('# P0 IOX ',(1x,e12.5,i5)) -906 format('# fact_B fact_T fact_n iscal',(3(1x,e12.5),i5)) -905 format('# EQUILIBRIUM CASE : ',a24) -907 format('# PROFILE file : ',a24) -909 format('# LAUNCHER CASE : ',a24) -910 format('# ipec nd ipsinorm sspl psdbnd : ',3i5,2(1x,e12.5)) -914 format('# |BT0| , sgnB_phi, sgnI_phi, icocos : ',3(1x,e12.5),i5) -915 format('# psi_1.5 rhop_1.5 rhot_1.5 : ',3(1x,e12.5)) -920 format('# psi_2 rhop_2 rhot_2 : ',3(1x,e12.5)) -911 format('# fghz : ',(1x,e12.5)) -912 format('# Cohen model') -913 format('# Momentum conservation model') +900 format('# Nray ktx rmx : ',2i5,1x,es12.5) +901 format('# igrad iwarm ilarm ieccd idst : ',6i5) +902 format('# X Y Z launching point (cm) : ',3(1x,es12.5)) +903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,es12.5)) +904 format('# GRAY revision : ',a) +905 format('# EQUILIBRIUM file : ',a24) +906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5)) +907 format('# PROFILES file : ',a24) +908 format('# alpha beta launching angles (deg) CYL : ',2(1x,es12.5)) +909 format('# LAUNCHER file : ',a24) +910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5) +911 format('# fghz P0 IOX : ',2(1x,es12.5),i5) +912 format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5)) +913 format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : ' + . ,8(1x,es12.5)) +914 format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5) +915 format('# dst nstep : ',1x,es12.5,i5) +916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2)) 99 format(20(1x,e16.8e3)) end @@ -871,7 +887,7 @@ c c c print circular magnetic surfaces iequil=1 c - write(71,*) '#i psi R z' + write(71,*) '#i psin R z' dal=2.0d-2*pi drn=0.1d0 do i=1,10 @@ -1305,7 +1321,7 @@ c call points_ox(rmaxis,zmaxis,rmop,zmop,psinop,info) rmaxis=rmop zmaxis=zmop - print'(a,2f8.4,e12.5)','O-point',rmop,zmop,psinop + print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinop c c search for X-point if ixp not = 0 c @@ -1315,7 +1331,7 @@ c z10=zbmin call points_ox(r10,z10,rxp,zxp,psinxp,info) if(psinxp.ne.-1.0d0) then - print'(a,2f8.4,e12.5)','X-point',rxp,zxp,psinxp + print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxp rbmin=rxp zbmin=zxp delpsinox=(psinxp-psinop) @@ -1481,7 +1497,6 @@ c tol = sqrt(dpmpar(1)) call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then -c print*,' subr points_tgo: info=',info end if rf=xvec(1) zf=xvec(2) @@ -1522,7 +1537,7 @@ c common/eqnn/nr,nz,npp,nintp common/dens/dens,ddens c - write(55,*) ' #psi rhop rhot ne Te q Jphi' + write(55,*) ' #psi rhot ne Te q Jphi' psin=0.0d0 rhop=0.0d0 rhot=0.0d0 @@ -1531,7 +1546,7 @@ c te=temp(psin) qq=qpsi(1) c - write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6 + write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 c nst=nr do i=2,nst @@ -1550,7 +1565,7 @@ c end if rhot=frhotor(psin) call tor_curr_psi(psin,ajphi) - write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6 + write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 end do c return @@ -2117,7 +2132,7 @@ c c computation of flux surface averaged quantities - write(71,*)' #i rhop R z' + write(71,*)' #i psin R z' nintp=nnintp ninpr=(nintp-1)/10 @@ -2333,7 +2348,7 @@ c end do end do - write(56,*)' #psi rhop rhot_eq rhot_av || |Bmx| |Bmn|'// + write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc' qqv(1)=qqv(2) @@ -2346,7 +2361,7 @@ c pstab(jp)=1.0d0 end if rhot_eq=frhotor_eq(pstab(jp)) - write(56,99) pstab(jp),rpstab(jp),rhot_eq,rhotqv(jp), + write(56,99) pstab(jp),rhot_eq,rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) end do @@ -2465,7 +2480,7 @@ c dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iov(jmx,kmx),tau1v(jmx,kmx) + dimension iiv(jmx,kmx),iov(jmx,kmx) common/iiv/iiv common/iov/iov @@ -2938,7 +2953,6 @@ c anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3) c if(abs(anpl).gt.0.99d0) then -c print*,' |Nparallel| > 1 !', sqrt(an2),anpl if(abs(anpl).le.1.05d0) then ierr=97 else @@ -3632,7 +3646,7 @@ c subroutine ic_gb implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36,zero=0.0d0,izero=0) + parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) @@ -3706,13 +3720,11 @@ c qi1=rci1-ui*ww1 qi2=rci2-ui*ww2 end if + c w01=sqrt(2.0d0*akinv/ww1) c z01=-rci1/(rci1**2+ww1**2) c w02=sqrt(2.0d0*akinv/ww2) c z02=-rci2/(rci2**2+ww2**2) -c -c spostare nel do???? -c qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2 @@ -3724,15 +3736,6 @@ c rciyy=dble(qqyy) rcixy=dble(qqxy)/2.0d0 -C d2ww1=-2.0d0*(dww1*rci1+ww1*drci1) -C d2ww2=-2.0d0*(dww2*rci2+ww2*drci2) -C d2rci1=2.0d0*(ww1*dww1-rci1*drci1) -C d2rci2=2.0d0*(ww2*dww2-rci2*drci2) -C dqi1=drci1-ui*dww1 -C dqi2=drci2-ui*dww2 -C d2qi1=d2rci1-ui*d2ww1 -C d2qi2=d2rci2-ui*d2ww2 -c dqi1=-qi1**2 dqi2=-qi2**2 d2qi1=2*qi1**3 @@ -3743,7 +3746,6 @@ c d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic) -c dwwxx=-dimag(dqqxx) dwwyy=-dimag(dqqyy) dwwxy=-dimag(dqqxy)/2.0d0 @@ -3782,7 +3784,6 @@ c dx0= x0t*csps+snps*(y0t*csth+z0t*snth) dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) dz0= z0t*csth-y0t*snth -c x0=x00+dx0 y0=y00+dy0 z0=z00+dz0 @@ -3885,7 +3886,13 @@ c write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . -1.0d0,zero,zero,zero end if - if(j.eq.1) write(17,99) zero,zero,zero,zero + if(j.eq.1.and.k.eq.1) then + write(17,99) zero,zero,zero,zero + write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + . zero,zero,zero,zero, + . zero,zero,zero,zero,one,zero,zero, + . zero,zero,one,zero,zero,zero,zero,one + end if end do end do c @@ -4034,7 +4041,13 @@ c write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . -1.0d0,zero,zero,zero end if - if(j.eq.1) write(17,99) zero,zero,zero,zero + if(j.eq.1.and.k.eq.1) then + write(17,99) zero,zero,zero,zero + write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + . zero,zero,zero,zero, + . zero,zero,zero,zero,one,zero,zero, + . zero,zero,one,zero,zero,zero,zero,one + end if end do end do c @@ -4053,7 +4066,7 @@ c - subroutine ic_rt0 + subroutine ic_rt2 implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36,zero=0.0d0,izero=0) @@ -4118,20 +4131,22 @@ c vgradi=0.0d0 ddi=2.0d0*vgradi c - if(j.eq.nray.or.j.eq.1) then + if(j.eq.nray) then r0=sqrt(x0**2+y0**2) x0m=x0/1.0d2 y0m=y0/1.0d2 r0m=r0/1.0d2 z0m=z0/1.0d2 - if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m - * ,zero,zero,zero,zero,zero - if(j.eq.1.and.k.eq.1) - * write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero - * ,zero,zero,zero + write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, + . -1.0d0,zero,zero,zero end if - if(k.eq.1.and.j.eq.nray) - * write(27,99) zero,x0,y0,z0,r0,dd,vgradi + if(j.eq.1.and.k.eq.1) then + write(17,99) zero,zero,zero,zero + write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + . zero,zero,zero,zero, + . zero,zero,zero,zero,one,zero,zero, + . zero,zero,one,zero,zero,zero,zero,one + end if end do end do c @@ -5627,8 +5642,7 @@ c end if if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti,yti,zti,rti, - . sqrt(psinv11) + . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 c if(rti.ge.rtimx.and.j.eq.nray) rtimx=rti if(rti.le.rtimn.and.j.eq.nray) rtimn=rti @@ -5636,14 +5650,13 @@ c end do c if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1, - . sqrt(psinv11) + . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 if(iproj.eq.1) write(nfile,*) ' ' end do c write(nfile,*) ' ' c - write(12,99) istep,st,sqrt(abs(psinv11)),rtimn,rtimx + write(12,99) istep,st,psinv11,rtimn,rtimx return 99 format(i5,12(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) @@ -5956,41 +5969,32 @@ c dPdV [MW/m^3], Jcd [MA/m^2] if(ieccd.eq.0) currt=0.0d0 currtka=currt*1.0d3 - if(index_rt.eq.1) then - write(7,*)'#B0 beta alpha Icd Pa Jphi '// - .'rhotj drhotj rhotjav rhotjava drhotjava dpdvmx rhotp drhotp '// - .'rhotpav drhotpav taumn taumx stmx polpsi polchi index_rt' - write(48,*) '#beta alpha rhop rhot dPdV Jphi Jcda P% Pins'// - . ' Icdins index_rt' - end if - write(6,*)' ' - write(6,*)'#beta alpha Icd Pa Jphi rhotj drhotj '// - .'rhotjav rhotjava drhotjava dpdvmx rhotp drhotp rhotpav '// - .'drhotpav taumn taumx stmx Pins_02 Pins_05 Pins_085' - write(6,99) betac,alfac,currtka,pabstot,ajmxfi,rhotjfi,drhotjfi, - . rhotjav,rhotjava,drhotjava, - . dpdvmx,rhotp,drhotp,rhotpav,drhotpav, - . taumn,taumx,stmx,pins_02,pins_05,pins_085 + write(6,*)'#beta alpha Icd Pa Jphimx dPdVmx '// + .'rhotj rhotjava rhotp rhotpav '// + .'drhotjava drhotpav stmx polpsi polchi index_rt' + write(6,99) betac,alfac,currtka,pabstot,ajmxfi,dpdvmx, + . rhotjfi,rhotjava,rhotp,rhotpav, + . drhotjava,drhotpav, + . stmx,psipol,chipol,real(index_rt) - write(7,99) btrcen,betac,alfac,currtka,pabstot, - . ajmxfi,rhotjfi, - . drhotjfi,rhotjav,rhotjava,drhotjava, - . dpdvmx,rhotp,drhotp,rhotpav,drhotpav, - . taumn,taumx,stmx,psipol,chipol,real(index_rt) + write(7,99) currtka,pabstot,ajmxfi,dpdvmx, + . rhotjfi,rhotjava,rhotp,rhotpav, + . drhotjava,drhotpav, + . stmx,psipol,chipol,real(index_rt) do i=1,nd if (ipec.eq.0) then - psip=rtab(i) + psin=rtab(i) rhop=sqrt(rtab(i)) else - psip=rtab(i)**2 + psin=rtab(i)**2 rhop=rtab(i) end if pinsr=0.0d0 if(pabstot.gt.0) pinsr=pins(i)/pabstot - write(48,99) betac,alfac,rhop,rhotv(i),dpdv(i),ajphiv(i) - . ,ajcdv(i),pinsr,pins(i),currins(i),real(index_rt) + write(48,99) psin,rhotv(i),ajphiv(i),ajcdv(i),dpdv(i), + . currins(i),pins(i),pinsr,real(index_rt) end do return @@ -6246,12 +6250,6 @@ c endif gam=atan(sngam/csgam)*180.d0/pi -c if(ipolc.eq.0.or.ipolc.eq.1) then -c write(28,111) beta,alfa,anpl,-gam,ffo,ffx,xe2om,xe2xm -c write(29,111) beta,alfa,qqom,uuom,vvom,psiom,chiom, -c . qqxm,uuxm,vvxm,psixm,chixm -c end if - return 111 format(20(1x,e12.5)) end