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

@@ -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)
+      subroutine eccd(yg,anpl,anprre,amu,zeff,
-      use green_func_p    
+     *  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)
+      real*8 cst2,eccdpar(5)
-      integer ieccd,nharm,nhf,nhn
+      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
+        do nhn=nhmn,nhmx
-          call curr_int(fjch,eccdpar,5,resj,resp)
+          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
+        do nhn=nhmn,nhmx
-          call curr_int(fjch0,eccdpar,1,resj,resp)
+          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
+        do nhn=nhmn,nhmx
-          call curr_int(fjncl,eccdpar,5,resj,resp)
+          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
+        end do
-      nhn=nhn-1
        
        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)