diff --git a/src/gray.f b/src/gray.f index e490d01..fd4a4f6 100644 --- a/src/gray.f +++ b/src/gray.f @@ -5939,13 +5939,11 @@ c apar(15) = ygn c npar=nfpp+necp - do i=1,necp - apar(nfpp+i) = eccdpar(i) - end do + apar(nfpp+1:npar) = eccdpar(1:necp) c if(duu.gt.1.0e-6_wp_) then - call dqagsmv(fpp,uu1,uu2,apar,npar,epsa,epsr,resp,epp,neval, - . ier,liw,lw,last,iw,w) + call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, + . epp,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) ierr=90 end if c @@ -6058,9 +6056,9 @@ c anpl=extrapar(2) amu=extrapar(3) anprre=extrapar(4) - ex=dcmplx(extrapar(5),extrapar(6)) - ey=dcmplx(extrapar(7),extrapar(8)) - ez=dcmplx(extrapar(9),extrapar(10)) + ex=cmplx(extrapar(5),extrapar(6),wp_) + ey=cmplx(extrapar(7),extrapar(8),wp_) + ez=cmplx(extrapar(9),extrapar(10),wp_) ithn=int(extrapar(11)) nhn=int(extrapar(12)) uplp=extrapar(13) @@ -6154,10 +6152,6 @@ c local variables c anpl=extrapar(2) anprre=extrapar(4) - ex=dcmplx(extrapar(5),extrapar(6)) - ey=dcmplx(extrapar(7),extrapar(8)) - ez=dcmplx(extrapar(9),extrapar(10)) - nhn=int(extrapar(12)) uplp=extrapar(13) uplm=extrapar(14) ygn=extrapar(15) @@ -6238,11 +6232,6 @@ c local variables complex(wp_) :: ex,ey,ez c anpl=extrapar(2) - anprre=extrapar(4) - ex=dcmplx(extrapar(5),extrapar(6)) - ey=dcmplx(extrapar(7),extrapar(8)) - ez=dcmplx(extrapar(9),extrapar(10)) - nhn=int(extrapar(12)) ygn=extrapar(15) zeff=extrapar(16) @@ -6291,28 +6280,35 @@ c extrapar(15) = ygn c c extrapar(16) = zeff c extrapar(17) = fc -c extrapar(18) = hb +c extrapar(18) = rbx c extrapar(19:18+(npar-18)/2) = tlm c extrapar(19+(npar-18)/2:npar) = chlm c use const_and_precisions, only : wp_ use green_func_p, only: GenSpitzFunc + use dierckx, only : splev,splder implicit none c arguments integer :: npar real(wp_) :: upl,fjncl real(wp_), dimension(npar) :: extrapar +c local constants + integer, parameter :: ksp=3 c local variables integer :: nlm - real(wp_) :: anpl,amu,ygn,zeff,fc,hb,gam,u2,u,upr2, + real(wp_) :: anpl,amu,ygn,zeff,fc,rbx,gam,u2,u,upr2, . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp +c local variables + integer :: ier + real(wp_), dimension((npar-18)/2) :: wrk + real(wp_), dimension(1) :: xs,ys c anpl=extrapar(2) amu=extrapar(3) ygn=extrapar(15) zeff=extrapar(16) fc=extrapar(17) - hb=extrapar(18) + rbx=extrapar(18) gam=anpl*upl+ygn u2=gam*gam-1.0_wp_ @@ -6324,13 +6320,18 @@ c fk=fk*(4.0_wp_/amu**2) dfk=dfk*(2.0_wp_/amu)*bth - alam=upr2/u2/hb + alam=upr2/u2/rbx nlm=(npar-18)/2 - call vlambda(alam,extrapar(19), - . extrapar(19+nlm), - . nlm,fu,dfu) + xs(1)=alam +c + call splev(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3, + . xs(1),ys(1),1,ier) + fu=ys(1) + call splder(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3,1, + . xs(1),ys(1),1,wrk,ier) + dfu=ys(1) - eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb + eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/rbx if(upl.lt.0) eta=-eta fjncl=eta*fpp(upl,extrapar,npar) return @@ -6338,32 +6339,32 @@ c c c c - subroutine vlambda(alam,tlm,chlm,nlmt,fv,dfv) - use const_and_precisions, only : wp_ - use dierckx, only : splev,splder - implicit none -c local constants - integer, parameter :: ksp=3 -c arguments - integer :: nlmt - real(wp_) :: alam,fv,dfv - real(wp_), dimension(nlmt) :: tlm,chlm -c local variables - integer :: nlm,ier - real(wp_), dimension(nlmt) :: wrk - real(wp_), dimension(1) :: xxs,ffs -c - nlm=nlmt - xxs(1)=alam -c - call splev(tlm,nlm,chlm,ksp,xxs(1),ffs(1),1,ier) - fv=ffs(1) -c - call splder(tlm,nlm,chlm,ksp,1,xxs(1),ffs(1),1,wrk,ier) - dfv=ffs(1) -c - return - end +C subroutine vlambda(alam,tlm,chlm,nlmt,fv,dfv) +C use const_and_precisions, only : wp_ +C use dierckx, only : splev,splder +C implicit none +Cc local constants +C integer, parameter :: ksp=3 +Cc arguments +C integer :: nlmt +C real(wp_) :: alam,fv,dfv +C real(wp_), dimension(nlmt) :: tlm,chlm +Cc local variables +C integer :: nlm,ier +C real(wp_), dimension(nlmt) :: wrk +C real(wp_), dimension(1) :: xxs,ffs +Cc +C nlm=nlmt +C xxs(1)=alam +Cc +C call splev(tlm,nlm,chlm,ksp,xxs(1),ffs(1),1,ier) +C fv=ffs(1) +Cc +C call splder(tlm,nlm,chlm,ksp,1,xxs(1),ffs(1),1,wrk,ier) +C dfv=ffs(1) +Cc +C return +C end c c c