diff --git a/src/gray.f b/src/gray.f index bc15cd2..b976a81 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1848,7 +1848,7 @@ c function frhotor_av(psi) implicit real*8(a-h,o-z) - parameter(nnintp=41) + parameter(nnintp=101) dimension rpstab(nnintp),crhotq(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,npp,nintp @@ -2117,21 +2117,6 @@ c 111 format(i6,12(1x,e12.5)) 99 format(12(1x,e12.5)) end - - subroutine sort(n,a) - implicit none - integer n,i,j - double precision a(n),temp - do i = 2, n - j = i - 1 - temp = a(i) - do while (j.ge.1.and.a(j).gt.temp) - a(j+1) = a(j) - j = j - 1 - end do - a(j+1) = temp - end do - end c c c @@ -2139,7 +2124,7 @@ c implicit real*8 (a-h,o-z) real*8 lam - parameter(nnintp=41,ncnt=100,ncntt=2*ncnt+1,nlam=41) + parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) parameter(zero=0.0d0,one=1.0d0) parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) @@ -2360,6 +2345,8 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qqv(jp)=qq +c +c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c computation of rhot from calculated q profile @@ -3963,12 +3950,15 @@ c . -1.0d0,zero,zero,zero,one end if 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 + ddr110=dd end if + if(j.eq.nrayr.and.k.eq.1) then + write(17,99) zero,ddr110,dd,ddi + end if end do end do c @@ -4282,7 +4272,7 @@ c subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, . bmxi,bmni,fci,intp) implicit real*8 (a-h,o-z) - parameter(nnintp=41) + parameter(nnintp=101) dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4) dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4) dimension carea(nnintp,4),cfc(nnintp,4) @@ -4316,7 +4306,7 @@ c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) implicit real*8 (a-h,o-z) - parameter(nnintp=41) + parameter(nnintp=101) dimension rpstab(nnintp) dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4) common/pstab/rpstab @@ -4660,6 +4650,7 @@ c common/resah/anpl2,dnl c common/cri/cr,ci + common/warm/iwarm,ilarm c anpl2=anpl**2 dnl=1.0d0-anpl2 @@ -4676,7 +4667,8 @@ c end do end do c - call hermitian(rr,lrm) + if(iwarm.eq.2) call hermitian(rr,lrm) + if(iwarm.eq.4) call hermitian_2(rr,lrm) c call antihermitian(ri,lrm) c @@ -4898,6 +4890,192 @@ c c return end + + + + subroutine hermitian_2(rr,lrm) + implicit real*8(a-h,o-z) + parameter(tmax=5,npts=500) + dimension rr(-lrm:lrm,0:2,0:lrm) + parameter(epsa=0.0d0,epsr=1.0d-4) + parameter (lw=5000,liw=lw/4) + dimension w(lw),iw(liw) + external fhermit +c + common/ygyg/yg + common/amut/amu + common/nplr/anpl,anpr + common/warm/iwarm,ilarm + common/cri/cr,ci + common/nmhermit/n,m,ih +c + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=0.0d0 + end do + end do + end do +c + llm=min(3,lrm) +c + bth2=2.0d0/amu + bth=sqrt(bth2) + amu2=amu*amu + amu4=amu2*amu2 + amu6=amu4*amu2 +c + n1=1 + if(iwarm.gt.10) n1=-llm +c + do n=n1,llm + nn=abs(n) + do m=nn,llm + if(n.eq.0.and.m.eq.0) then + ih=2 +c call dqagi(fhermit,bound,2,epsa,epsr,resfh, + call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, + . epp,neval,ier,liw,lw,last,iw,w) + rr(n,2,m) = resfh + else + do ih=0,2 +c call dqagi(fhermit,bound,2,epsa,epsr,resfh, + call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, + . epp,neval,ier,liw,lw,last,iw,w) + rr(n,ih,m) = resfh + end do + end if + end do + end do +c write(83,'(12(1x,e12.5))') +c . yg,rr(1,0,1),rr(1,1,1),rr(1,2,1) + + if(iwarm.gt.10) return +c + sy1=1.0d0+yg + sy2=1.0d0+yg*2.0d0 + sy3=1.0d0+yg*3.0d0 +c + bth4=bth2*bth2 + bth6=bth4*bth2 +c + anpl2=anpl*anpl +c + rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) + . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) + rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) + rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) + rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 + . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ + . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) + rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ + . 1.5d0*anpl2/sy1**2)) + rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* + . anpl2/sy1**2)) +c + if(llm.gt.1) then +c + rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ + . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) + rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) + rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) + rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* + . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* + . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 + . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) + rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* + . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) + rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* + . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) + rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* + . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* + . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 + . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) + rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* + . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) + rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* + . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) +c + if(llm.gt.2) then +c + rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* + . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) + rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) + rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) + rr(-1,0,3) = -12.0d0*bth4/sy1* + . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ + . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) + . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) + rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* + . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) + rr(-1,2,3) = -6.0d0*bth6/sy1* + . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) +c + rr(-2,0,3) = -12.0d0*bth4/sy2* + . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ + . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) + . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) + rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* + . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) + rr(-2,2,3) = -6.0d0*bth6/sy2* + . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) +c + rr(-3,0,3) = -12.0d0*bth4/sy3* + . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ + . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) + . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) + rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* + . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) + rr(-3,2,3) = -6.0d0*bth6/sy3* + . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) +c + end if + end if +c + return + end + + + + function fhermit(t) + implicit real*8 (a-h,o-z) +c + common/ygyg/yg + common/amut/amu + common/nplr/anpl,anpr + common/cri/cr,ci + common/nmhermit/n,m,ih + + bth2=2.0d0/amu + bth=sqrt(bth2) + amu2=amu*amu + amu4=amu2*amu2 + amu6=amu4*amu2 + + rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=1.0d0+t*t/amu + exdxdt=cr*exp(-t*t)*gx/rxt + nn=abs(n) + gr=anpl*upl+n*yg + zm=-amu*(gx-gr) + s=amu*(gx+gr) + zm2=zm*zm + zm3=zm2*zm + call calcei3(zm,fe0m) + ffe=0.0d0 + uplh=upl**ih + if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 + if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2 + if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s + . +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4 + if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) + . +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6 + fhermit= exdxdt*ffe + return + end c c c @@ -5626,7 +5804,7 @@ c gg=F(u)/u with F(u) as in Cohen paper subroutine vlambda(alam,psi,fv,dfv) implicit real*8 (a-h,o-z) - parameter (nnintp=41,nlam=41) + parameter (nnintp=101,nlam=41) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54)