istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

few fixes in pec plus other minor changes

Browse Source
This commit is contained in:
Lorenzo Figini 2015-02-18 16:17:21 +00:00
parent 608d63acfe
commit 1df304f221
1 changed files with 36 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
src/gray.f
View File

@ -69,7 +69,12 @@ c second pass into plasma
subroutine gray_integration subroutine gray_integration
implicit real*8 (a-h,o-z) implicit none
integer istep,istop,index_rt
integer nstep
real*8 st,dst,strfl11
integer i
real*8 st0
common/ss/st common/ss/st
common/dsds/dst common/dsds/dst
@ -137,24 +142,24 @@ c
c c
c print final results on screen c print final results on screen
c c
print*,' ' write(*,*)
print'(a,f9.4)','final step (s, ct, Sr) = ',st write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st
if(iwarm.gt.0) then if(iwarm.gt.0) then
print '(a,2e12.5)','taumn, taumx = ', taumn,taumx write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx
else else
print '(a,2f9.4)','taumn, taumx = ', zero,zero write(*,'(a,2f9.4)') 'taumn, taumx = ', zero,zero
end if end if
c c
print'(a,f9.4)','Pabs_tot (MW) = ',pabstot write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot
currtka =currtot*1.0d3 currtka =currtot*1.0d3
print'(a,f9.4)','I_tot (kA) = ',currtka write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka
c c
if (index_rt.eq.1) then if (index_rt.eq.1) then
if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq
if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM' if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM'
if(iprof.eq.1) write(6,*) 'PROFILE file : ',filenmprf if(iprof.eq.1) write(*,*) 'PROFILE file : ',filenmprf
if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0
if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm if(ibeam.ge.1) write(*,*) 'LAUNCHER CASE : ',filenmbm
end if end if
c c
@ -248,7 +253,7 @@ c
call gwork(j,k) call gwork(j,k)
c c
if(ierr.gt.0) then if(ierr.gt.0) then
print*,' IERR = ', ierr write(*,*) ' IERR = ', ierr
if(ierr.eq.97) then if(ierr.eq.97) then
c igrad=0 c igrad=0
c ierr=0 c ierr=0
@ -959,9 +964,9 @@ c set simple limiter as two cylindrical walls at rwallm and r00
rlim(3)=rlim(2) rlim(3)=rlim(2)
rlim(4)=rlim(1) rlim(4)=rlim(1)
rlim(5)=rlim(1) rlim(5)=rlim(1)
zlim(1)=-dst*nmx*1.d-2 zlim(1)=(z00-dst*nmx)*1.d-2
zlim(2)=zlim(1) zlim(2)=zlim(1)
zlim(3)=dst*nmx*1.d-2 zlim(3)=(z00+dst*nmx)*1.d-2
zlim(4)=zlim(3) zlim(4)=zlim(3)
zlim(5)=zlim(1) zlim(5)=zlim(1)
ipass=abs(ipass) ipass=abs(ipass)
@ -1568,8 +1573,8 @@ c compute B_toroidal on axis
btaxis=fpol(1)/rmaxis btaxis=fpol(1)/rmaxis
btrcen=abs(btrcen)*factb btrcen=abs(btrcen)*factb
print'(a,f8.4)','factb = ',factb print'(a,f8.4)','factb = ',factb
print'(a,f8.4)','|BT_centr|= ',btrcen print'(a,f8.4)','BT_centr= ',btrcen
print'(a,f8.4)','|BT_axis| = ',abs(btaxis) print'(a,f8.4)','BT_axis = ',btaxis
c compute normalized rho_tor from eqdsk q profile c compute normalized rho_tor from eqdsk q profile
call rhotor(nr) call rhotor(nr)
@ -3693,8 +3698,8 @@ c
call fpbisp(tr,nsr,tz,nsz,ccspl,3,3, call fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
. rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
psinv=(ffspl(1)-psinop)/psiant psinv=(ffspl(1)-psinop)/psiant
c if(psinv.lt.0.0d0) if(psinv.lt.0.0d0)
c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim
c c
nur=1 nur=1
nuz=0 nuz=0
@ -3945,8 +3950,7 @@ c
fzeff=zeff fzeff=zeff
else else
call locate(psrad,npp,ps,k) call locate(psrad,npp,ps,k)
if(k.eq.0) k=1 k=max(1,min(k,npp-1))
if(k.eq.npp) k=npp-1
dps=ps-psrad(k) dps=ps-psrad(k)
fzeff=spli(cz,npmx,k,dps) fzeff=spli(cz,npmx,k,dps)
endif endif
@ -6420,7 +6424,8 @@ c
else else
if(xxi(i).gt.rtbc) exit if(xxi(i).gt.rtbc) exit
if(xxi(i).lt.xxi(i-1)) then if(xxi(i).lt.xxi(i-1)) then
isev(is)=i ! isev(is)=i !!!!!!!!!! it should be isev(is)=i-1
isev(is)=i-1
is=is+1 is=is+1
idecr=-1 idecr=-1
end if end if
@ -6446,12 +6451,13 @@ c
ind2=nd ind2=nd
iind=1 iind=1
end if end if
do ind=ind1,ind2,iind do ind=ind1,ind2,iind !!!!!!!!!! is it safe?
iind=ind ! iind=ind !!!!!!!!!! iind reused in the loop!
if (idecr.eq.-1) iind=ind-1 indi=ind !!!!!!!!!! iind reused in the loop!
rt1=rtab1(iind) if (idecr.eq.-1) indi=ind-1
rt1=rtab1(indi)
call locatex(xxi,iise,iise0,iise,rt1,itb1) call locatex(xxi,iise,iise0,iise,rt1,itb1)
if(itb1.gt.iise0.and.itb1.lt.iise) then if(itb1.ge.iise0.and.itb1.lt.iise) then
call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1), call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),
. ypt(itb1+1),rt1,ppa2) . ypt(itb1+1),rt1,ppa2)
call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1), call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),
@ -6534,12 +6540,12 @@ c of gaussian profile
rhpp=frhopol(rhotpav) rhpp=frhopol(rhotpav)
rhpj=frhopol(rhotjava) rhpj=frhopol(rhotjava)
dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp))
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
. drhotp,drhopp)
if(ieccd.ne.0) then
ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
. drhotp,drhopp)
if(ieccd.ne.0) then
call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
. drhotjfi,drhopfi) . drhotjfi,drhopfi)
xps=rhopfi xps=rhopfi