From 045c581865c70bee90ff4d1b1c619851f559ee03 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Wed, 13 May 2015 16:27:52 +0000 Subject: [PATCH] other variables moved from common to argument list in eccd. fixed print of nharm in print_output, added missin /ierr/ common in ecrh_cd --- src/gray.f | 91 ++++++++++++++++++++++++------------------------------ 1 file changed, 40 insertions(+), 51 deletions(-) diff --git a/src/gray.f b/src/gray.f index 6bd1aa9..86c18cc 100644 --- a/src/gray.f +++ b/src/gray.f @@ -489,7 +489,7 @@ c common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst - common/nhn/nhn + common/nharm/nharm,nhf common/iokh/iohkawa common/p0/p0mw common/tau1v/tau1v @@ -584,7 +584,7 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 . psi,rhot,dens11,tekev, . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, . alpha11,tau11,pt11,ajcd11,dids11, - . dble(nhn),dble(iohkawa),dble(index_rt) + . dble(nhf),dble(iohkawa),dble(index_rt) c call polarcold(exf,eyif,ezf,elf,etf) end if @@ -4771,7 +4771,6 @@ c parameter(vcsi=2.99792458d+8) parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3) c - common/ithn/ithn common/nharm/nharm,nhf common/warm/iwarm,ilarm common/ieccd/ieccd @@ -4784,6 +4783,9 @@ c common/nprw/anprre,anprim c common/absor/alpha,effjcd,akim,tau +c + common/ierr/ierr + common/iokh/iokhawa c common/psival/psinv common/tete/tekev @@ -4850,7 +4852,8 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm c ithn=1 if(lrm.gt.nharm) ithn=2 - if(ieccd.gt.0) call eccd(effjcd) + if(ieccd.gt.0) call eccd(yg,anpl,anprre,amu,zeff, + * ieccd,nharm,nhf,ithn,effjcd,iokhawa,ierr) return 999 format(12(1x,e12.5)) @@ -5793,45 +5796,45 @@ c - subroutine eccd(effjcd) - use green_func_p + subroutine eccd(yg,anpl,anprre,amu,zeff, + * ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) + use green_func_p, only: SpitzFuncCoeff,pi implicit none real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev - real*8 anucc,canucc,ddens + real*8 yg,anpl,anpr,amu,anprre,anprim + real*8 mc2,anucc,canucc,ddens real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff,psinv real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc real*8 fjch,fjncl,fjch0,fconic - real*8 cst,cst2,eccdpar(5) - integer ieccd,nharm,nhf,nhn + real*8 cst2,eccdpar(5) + complex*16 ex,ey,ez + integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa external fjch,fjncl,fjch0 parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(vcsi=2.99792458d+8) parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2) - parameter(canucc=4.0d13*pi*qe**4/(me**2*vc**3)) + parameter(mc2=mesi*vcsi*vcsi/qesi*1d-3) + parameter(canucc=2.0d13*pi*qe**4/(me**2*vc**3)) parameter(ceff=qesi/(mesi*vcsi)) - common/nharm/nharm,nhf - common/nhn/nhn - - common/ieccd/ieccd - common/tete/tekev - common/dens/dens,ddens - common/zz/Zeff common/psival/psinv + + common/dens/dens,ddens + common/epolar/ex,ey,ez common/btot/btot common/bmxmn/bmax,bmin common/fc/fc - common/cst/cst,cst2 anum=0.0d0 denom=0.0d0 effjcd=0.0d0 - coullog=24.0d0-log(1.0d4*sqrt(0.1d0*dens)/tekev) + tekev=mc2/amu + coullog=48.0d0-log(1.0d7*dens/tekev**2) anucc=canucc*dens*coullog -c nhf=nharm +c nhmx=nhmn eccdpar(1)=zeff @@ -5849,7 +5852,6 @@ c fp0s= P_a (alams) rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 - cst=sqrt(cst2) alams=sqrt(1.0d0-bmin/bmax) pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0 fp0s=fconic(alams,pa,0) @@ -5857,17 +5859,18 @@ c fp0s= P_a (alams) eccdpar(3)=alams eccdpar(4)=pa eccdpar(5)=fp0s - do nhn=nharm,nhf - call curr_int(fjch,eccdpar,5,resj,resp) + do nhn=nhmn,nhmx + call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, + * fjch,eccdpar,5,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do case(2:9) - cst=0.0d0 cst2=0.0d0 - do nhn=nharm,nhf - call curr_int(fjch0,eccdpar,1,resj,resp) + do nhn=nhmn,nhmx + call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, + * fjch0,eccdpar,1,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do @@ -5880,18 +5883,17 @@ c rzfc=(1+Zeff)/fc rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 - cst=sqrt(cst2) call SpitzFuncCoeff(Tekev,Zeff,fc) eccdpar(2) = tekev eccdpar(3) = fc eccdpar(4) = rbx eccdpar(5) = psinv - do nhn=nharm,nhf - call curr_int(fjncl,eccdpar,5,resj,resp) + do nhn=nhmn,nhmx + call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, + * fjncl,eccdpar,5,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp - end do - nhn=nhn-1 + end do CASE DEFAULT print*,'ieccd undefined' @@ -5906,30 +5908,17 @@ c c c c - subroutine curr_int(fcur,eccdpar,necp,resj,resp) + subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, + * fcur,eccdpar,necp,resj,resp,iokhawa,ierr) implicit real*8(a-h,o-z) parameter(epsa=0.0d0,epsr=1.0d-2) parameter(xxcr=16.0d0) parameter (lw=5000,liw=lw/4) - parameter (npar=20) + parameter (nfpp=15) complex*16 ex,ey,ez - dimension w(lw),iw(liw),eccdpar(necp),apar(npar) + dimension w(lw),iw(liw),eccdpar(necp),apar(nfpp+necp) external fcur,fpp - common/ierr/ierr - common/iokh/iokhawa - common/cst/cst,cst2 - -c used and passed - common/ygyg/yg - common/nplr/anpl,anpr - common/amut/amu - common/nhn/nhn -c passed unused - common/epolar/ex,ey,ez - common/nprw/anprre,anprim - common/ithn/ithn - c EC power and current densities anpl2=anpl*anpl @@ -5978,8 +5967,9 @@ c apar(14) = uplm apar(15) = ygn c + npar=nfpp+necp do i=1,necp - apar(15+i) = eccdpar(i) + apar(nfpp+i) = eccdpar(i) end do c if(duu.gt.1.d-6) then @@ -6006,6 +5996,7 @@ c resonance curve crosses the trapping region c iokhawa=1 rdut=sqrt(rdu2t) + cst=sqrt(cst2) upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2) upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2) c @@ -6177,7 +6168,6 @@ c ex=dcmplx(extrapar(5),extrapar(6)) ey=dcmplx(extrapar(7),extrapar(8)) ez=dcmplx(extrapar(9),extrapar(10)) - ithn=int(extrapar(11)) nhn=int(extrapar(12)) uplp=extrapar(13) uplm=extrapar(14) @@ -6255,7 +6245,6 @@ c ex=dcmplx(extrapar(5),extrapar(6)) ey=dcmplx(extrapar(7),extrapar(8)) ez=dcmplx(extrapar(9),extrapar(10)) - ithn=int(extrapar(11)) nhn=int(extrapar(12)) ygn=extrapar(15) zeff=extrapar(16)