diff --git a/src/conical.f90 b/src/conical.f90 index ff7195b..fc9c606 100644 --- a/src/conical.f90 +++ b/src/conical.f90 @@ -57,7 +57,7 @@ contains return end if else - ti=cmplx(0._wp_,tau) + ti=cmplx(0._wp_,tau,wp_) ! if((-1._wp_ < x .and. x <= 0.0_wp_).or. & (0.0_wp_ < x .and. x <= 0.1_wp_ .and.tau<= 17.0_wp_).or. & @@ -316,7 +316,7 @@ contains return end if f=abs(t) - v=cmplx(x,f) + v=cmplx(x,f,wp_) if(x < 0.0_wp_) v=1.0_wp_-v h=(0.0_wp_,0.0_wp_) c=real(v) @@ -327,7 +327,7 @@ contains a=atan2(d,c) do i = 1,n c=c+1.0_wp_ - v=cmplx(c,d) + v=cmplx(c,d,wp_) h=h*v a=a+atan2(d,c) end do @@ -347,7 +347,7 @@ contains f=sin(c) e=d+0.5_wp_*log(e*f**2+0.25_wp_*(1.0_wp_-e)**2) f=atan2(cos(c)*tanh(d),f)-a*pi - clogam=1.1447298858494_wp_-cmplx(e,f)-clogam + clogam=1.1447298858494_wp_-cmplx(e,f,wp_)-clogam ! end if if(t < 0.0_wp_) clogam=conjg(clogam) diff --git a/src/eierf.f90 b/src/eierf.f90 index ef3dd8f..fe2fb58 100644 --- a/src/eierf.f90 +++ b/src/eierf.f90 @@ -629,278 +629,278 @@ contains result = ei end subroutine calcei3 - subroutine calerf(arg,result,jintt) -!------------------------------------------------------------------ +! subroutine calerf(arg,result,jintt) +!!------------------------------------------------------------------ +!! +!! this packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) +!! for a real argument x. it contains three function type +!! subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx), +!! and one subroutine type subprogram, calerf. the calling +!! statements for the primary entries are: +!! +!! y=erf(x) (or y=derf(x)), +!! +!! y=erfc(x) (or y=derfc(x)), +!! and +!! y=erfcx(x) (or y=derfcx(x)). +!! +!! the routine calerf is intended for internal packet use only, +!! all computations within the packet being concentrated in this +!! routine. the function subprograms invoke calerf with the +!! statement +!! +!! call calerf(arg,result,jintt) +!! +!! where the parameter usage is as follows +!! +!! function parameters for calerf +!! call arg result jintt +!! +!! erf(arg) any real argument erf(arg) 0 +!! erfc(arg) abs(arg) < xbig erfc(arg) 1 +!! erfcx(arg) xneg < arg < xmax erfcx(arg) 2 +!! +!!******************************************************************* +!!******************************************************************* +!! +!! Explanation of machine-dependent constants +!! +!! XMIN = the smallest positive floating-point number. +!! XINF = the largest positive finite floating-point number. +!! XNEG = the largest negative argument acceptable to ERFCX; +!! the negative of the solution to the equation +!! 2*exp(x*x) = XINF. +!! XSMALL = argument below which erf(x) may be represented by +!! 2*x/sqrt(pi) and above which x*x will not underflow. +!! A conservative value is the largest machine number X +!! such that 1.0 + X = 1.0 to machine precision. +!! XBIG = largest argument acceptable to ERFC; solution to +!! the equation: W(x) * (1-0.5/x**2) = XMIN, where +!! W(x) = exp(-x*x)/[x*sqrt(pi)]. +!! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to +!! machine precision. A conservative value is +!! 1/[2*sqrt(XSMALL)] +!! XMAX = largest acceptable argument to ERFCX; the minimum +!! of XINF and 1/[sqrt(pi)*XMIN]. +!! +!!******************************************************************* +!!******************************************************************* +!! +!! error returns +!! +!! the program returns erfc = 0 for arg >= xbig; +!! +!! erfcx = xinf for arg < xneg; +!! and +!! erfcx = 0 for arg >= xmax. +!! +!! +!! intrinsic functions required are: +!! +!! abs, aint, exp +!! +!! +!! author: w. j. cody +!! mathematics and computer science division +!! argonne national laboratory +!! argonne, il 60439 +!! +!! latest modification: march 19, 1990 +!! +!!------------------------------------------------------------------ +! implicit none +! real(wp_), intent(in) :: arg +! real(wp_), intent(out) :: result +! integer, intent(in) :: jintt +! integer :: i +! real(wp_) :: del,x,xden,xnum,y,ysq +!!------------------------------------------------------------------ +!! mathematical constants +!!------------------------------------------------------------------ +! real(wp_), parameter :: sqrpi=5.6418958354775628695e-1_wp_, & +! thresh=0.46875_wp_ +!!------------------------------------------------------------------ +!! machine-dependent constants +!!------------------------------------------------------------------ +! real(wp_), parameter :: xinf=1.79e308_wp_, & ! ~huge +! xneg=-26.628_wp_, & ! ? +! xsmall=1.11e-16_wp_, & ! ~epsilon/2 +! xbig=26.543_wp_, & ! ? +! xhuge=6.71e7_wp_, & ! ~1/sqrt(epsilon) +! xmax=2.53e307_wp_ ! ? +!!------------------------------------------------------------------ +!! coefficients for approximation to erf in first interval +!!------------------------------------------------------------------ +! real(wp_), dimension(5), parameter :: & +! a=(/3.16112374387056560_wp_,1.13864154151050156e02_wp_, & +! 3.77485237685302021e02_wp_,3.20937758913846947e03_wp_, & +! 1.85777706184603153e-1_wp_/) +! real(wp_), dimension(4), parameter :: & +! b=(/2.36012909523441209e01_wp_,2.44024637934444173e02_wp_, & +! 1.28261652607737228e03_wp_,2.84423683343917062e03_wp_/) +!!------------------------------------------------------------------ +!! coefficients for approximation to erfc in second interval +!!------------------------------------------------------------------ +! real(wp_), dimension(9), parameter :: & +! c=(/5.64188496988670089e-1_wp_,8.88314979438837594_wp_, & +! 6.61191906371416295e01_wp_,2.98635138197400131e02_wp_, & +! 8.81952221241769090e02_wp_,1.71204761263407058e03_wp_, & +! 2.05107837782607147e03_wp_,1.23033935479799725e03_wp_, & +! 2.15311535474403846e-8_wp_/) +! real(wp_), dimension(8), parameter :: & +! d=(/1.57449261107098347e01_wp_,1.17693950891312499e02_wp_, & +! 5.37181101862009858e02_wp_,1.62138957456669019e03_wp_, & +! 3.29079923573345963e03_wp_,4.36261909014324716e03_wp_, & +! 3.43936767414372164e03_wp_,1.23033935480374942e03_wp_/) +!!------------------------------------------------------------------ +!! coefficients for approximation to erfc in third interval +!!------------------------------------------------------------------ +! real(wp_), dimension(6), parameter :: & +! p=(/3.05326634961232344e-1_wp_,3.60344899949804439e-1_wp_, & +! 1.25781726111229246e-1_wp_,1.60837851487422766e-2_wp_, & +! 6.58749161529837803e-4_wp_,1.63153871373020978e-2_wp_/) +! real(wp_), dimension(5), parameter :: & +! q=(/2.56852019228982242_wp_,1.87295284992346047_wp_, & +! 5.27905102951428412e-1_wp_,6.05183413124413191e-2_wp_, & +! 2.33520497626869185e-3_wp_/) +!!------------------------------------------------------------------ +! x = arg +! y = abs(x) +! if (y <= thresh) then +!!------------------------------------------------------------------ +!! evaluate erf for |x| <= 0.46875 +!!------------------------------------------------------------------ +! ysq = zero +! if (y > xsmall) ysq = y * y +! xnum = a(5)*ysq +! xden = ysq +! do i = 1, 3 +! xnum = (xnum + a(i)) * ysq +! xden = (xden + b(i)) * ysq +! end do +! result = x * (xnum + a(4)) / (xden + b(4)) +! if (jintt /= 0) result = one - result +! if (jintt == 2) result = exp(ysq) * result +! return +!!------------------------------------------------------------------ +!! evaluate erfc for 0.46875 <= |x| <= 4.0 +!!------------------------------------------------------------------ +! else if (y <= four) then +! xnum = c(9)*y +! xden = y +! do i = 1, 7 +! xnum = (xnum + c(i)) * y +! xden = (xden + d(i)) * y +! end do +! result = (xnum + c(8)) / (xden + d(8)) +! if (jintt /= 2) then +! ysq = aint(y*sixten)/sixten +! del = (y-ysq)*(y+ysq) +! result = exp(-ysq*ysq) * exp(-del) * result +! end if +!!------------------------------------------------------------------ +!! evaluate erfc for |x| > 4.0 +!!------------------------------------------------------------------ +! else if (y < xbig .or. (y < xmax .and. jintt == 2)) then +! ysq = one / (y * y) +! xnum = p(6)*ysq +! xden = ysq +! do i = 1, 4 +! xnum = (xnum + p(i)) * ysq +! xden = (xden + q(i)) * ysq +! end do +! result = ysq *(xnum + p(5)) / (xden + q(5)) +! result = (sqrpi - result) / y +! if (jintt /= 2) then +! ysq = aint(y*sixten)/sixten +! del = (y-ysq)*(y+ysq) +! result = exp(-ysq*ysq) * exp(-del) * result +! end if +! else if (y >= xhuge) then +! result = sqrpi / y +! else +! result = zero +! end if +!!------------------------------------------------------------------ +!! fix up for negative argument, erf, etc. +!!------------------------------------------------------------------ +! if (jintt == 0) then +! result = (half - result) + half +! if (x < zero) result = -result +! else if (jintt == 1) then +! if (x < zero) result = two - result +! else +! if (x < zero) then +! if (x < xneg) then +! result = xinf +! else +! ysq = aint(x*sixten)/sixten +! del = (x-ysq)*(x+ysq) +! y = exp(ysq*ysq) * exp(del) +! result = (y+y) - result +! end if +! end if +! end if +! end subroutine calerf ! -! this packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) -! for a real argument x. it contains three function type -! subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx), -! and one subroutine type subprogram, calerf. the calling -! statements for the primary entries are: +! function derf(x) +!!-------------------------------------------------------------------- +!! +!! this subprogram computes approximate values for erf(x). +!! (see comments heading calerf). +!! +!! author/date: w. j. cody, january 8, 1985 +!! +!!-------------------------------------------------------------------- +! implicit none +! real(wp_) :: derf +! real(wp_), intent(in) :: x +! integer :: jintt +! real(wp_) :: result +!!------------------------------------------------------------------ +! jintt = 0 +! call calerf(x,result,jintt) +! derf = result +! end function derf ! -! y=erf(x) (or y=derf(x)), +! function derfc(x) +!!-------------------------------------------------------------------- +!! +!! this subprogram computes approximate values for erfc(x). +!! (see comments heading calerf). +!! +!! author/date: w. j. cody, january 8, 1985 +!! +!!-------------------------------------------------------------------- +! implicit none +! real(wp_) :: derfc +! real(wp_), intent(in) :: x +! integer :: jintt +! real(wp_) :: result +!!------------------------------------------------------------------ +! jintt = 1 +! call calerf(x,result,jintt) +! derfc = result +! end function derfc ! -! y=erfc(x) (or y=derfc(x)), -! and -! y=erfcx(x) (or y=derfcx(x)). -! -! the routine calerf is intended for internal packet use only, -! all computations within the packet being concentrated in this -! routine. the function subprograms invoke calerf with the -! statement -! -! call calerf(arg,result,jintt) -! -! where the parameter usage is as follows -! -! function parameters for calerf -! call arg result jintt -! -! erf(arg) any real argument erf(arg) 0 -! erfc(arg) abs(arg) < xbig erfc(arg) 1 -! erfcx(arg) xneg < arg < xmax erfcx(arg) 2 -! -!******************************************************************* -!******************************************************************* -! -! Explanation of machine-dependent constants -! -! XMIN = the smallest positive floating-point number. -! XINF = the largest positive finite floating-point number. -! XNEG = the largest negative argument acceptable to ERFCX; -! the negative of the solution to the equation -! 2*exp(x*x) = XINF. -! XSMALL = argument below which erf(x) may be represented by -! 2*x/sqrt(pi) and above which x*x will not underflow. -! A conservative value is the largest machine number X -! such that 1.0 + X = 1.0 to machine precision. -! XBIG = largest argument acceptable to ERFC; solution to -! the equation: W(x) * (1-0.5/x**2) = XMIN, where -! W(x) = exp(-x*x)/[x*sqrt(pi)]. -! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to -! machine precision. A conservative value is -! 1/[2*sqrt(XSMALL)] -! XMAX = largest acceptable argument to ERFCX; the minimum -! of XINF and 1/[sqrt(pi)*XMIN]. -! -!******************************************************************* -!******************************************************************* -! -! error returns -! -! the program returns erfc = 0 for arg >= xbig; -! -! erfcx = xinf for arg < xneg; -! and -! erfcx = 0 for arg >= xmax. -! -! -! intrinsic functions required are: -! -! abs, aint, exp -! -! -! author: w. j. cody -! mathematics and computer science division -! argonne national laboratory -! argonne, il 60439 -! -! latest modification: march 19, 1990 -! -!------------------------------------------------------------------ - implicit none - real(wp_), intent(in) :: arg - real(wp_), intent(out) :: result - integer, intent(in) :: jintt - integer :: i - real(wp_) :: del,x,xden,xnum,y,ysq -!------------------------------------------------------------------ -! mathematical constants -!------------------------------------------------------------------ - real(wp_), parameter :: sqrpi=5.6418958354775628695e-1_wp_, & - thresh=0.46875_wp_ -!------------------------------------------------------------------ -! machine-dependent constants -!------------------------------------------------------------------ - real(wp_), parameter :: xinf=1.79e308_wp_, & ! ~huge - xneg=-26.628_wp_, & ! ? - xsmall=1.11e-16_wp_, & ! ~epsilon/2 - xbig=26.543_wp_, & ! ? - xhuge=6.71e7_wp_, & ! ~1/sqrt(epsilon) - xmax=2.53e307_wp_ ! ? -!------------------------------------------------------------------ -! coefficients for approximation to erf in first interval -!------------------------------------------------------------------ - real(wp_), dimension(5), parameter :: & - a=(/3.16112374387056560_wp_,1.13864154151050156e02_wp_, & - 3.77485237685302021e02_wp_,3.20937758913846947e03_wp_, & - 1.85777706184603153e-1_wp_/) - real(wp_), dimension(4), parameter :: & - b=(/2.36012909523441209e01_wp_,2.44024637934444173e02_wp_, & - 1.28261652607737228e03_wp_,2.84423683343917062e03_wp_/) -!------------------------------------------------------------------ -! coefficients for approximation to erfc in second interval -!------------------------------------------------------------------ - real(wp_), dimension(9), parameter :: & - c=(/5.64188496988670089e-1_wp_,8.88314979438837594_wp_, & - 6.61191906371416295e01_wp_,2.98635138197400131e02_wp_, & - 8.81952221241769090e02_wp_,1.71204761263407058e03_wp_, & - 2.05107837782607147e03_wp_,1.23033935479799725e03_wp_, & - 2.15311535474403846e-8_wp_/) - real(wp_), dimension(8), parameter :: & - d=(/1.57449261107098347e01_wp_,1.17693950891312499e02_wp_, & - 5.37181101862009858e02_wp_,1.62138957456669019e03_wp_, & - 3.29079923573345963e03_wp_,4.36261909014324716e03_wp_, & - 3.43936767414372164e03_wp_,1.23033935480374942e03_wp_/) -!------------------------------------------------------------------ -! coefficients for approximation to erfc in third interval -!------------------------------------------------------------------ - real(wp_), dimension(6), parameter :: & - p=(/3.05326634961232344e-1_wp_,3.60344899949804439e-1_wp_, & - 1.25781726111229246e-1_wp_,1.60837851487422766e-2_wp_, & - 6.58749161529837803e-4_wp_,1.63153871373020978e-2_wp_/) - real(wp_), dimension(5), parameter :: & - q=(/2.56852019228982242_wp_,1.87295284992346047_wp_, & - 5.27905102951428412e-1_wp_,6.05183413124413191e-2_wp_, & - 2.33520497626869185e-3_wp_/) -!------------------------------------------------------------------ - x = arg - y = abs(x) - if (y <= thresh) then -!------------------------------------------------------------------ -! evaluate erf for |x| <= 0.46875 -!------------------------------------------------------------------ - ysq = zero - if (y > xsmall) ysq = y * y - xnum = a(5)*ysq - xden = ysq - do i = 1, 3 - xnum = (xnum + a(i)) * ysq - xden = (xden + b(i)) * ysq - end do - result = x * (xnum + a(4)) / (xden + b(4)) - if (jintt /= 0) result = one - result - if (jintt == 2) result = exp(ysq) * result - return -!------------------------------------------------------------------ -! evaluate erfc for 0.46875 <= |x| <= 4.0 -!------------------------------------------------------------------ - else if (y <= four) then - xnum = c(9)*y - xden = y - do i = 1, 7 - xnum = (xnum + c(i)) * y - xden = (xden + d(i)) * y - end do - result = (xnum + c(8)) / (xden + d(8)) - if (jintt /= 2) then - ysq = aint(y*sixten)/sixten - del = (y-ysq)*(y+ysq) - result = exp(-ysq*ysq) * exp(-del) * result - end if -!------------------------------------------------------------------ -! evaluate erfc for |x| > 4.0 -!------------------------------------------------------------------ - else if (y < xbig .or. (y < xmax .and. jintt == 2)) then - ysq = one / (y * y) - xnum = p(6)*ysq - xden = ysq - do i = 1, 4 - xnum = (xnum + p(i)) * ysq - xden = (xden + q(i)) * ysq - end do - result = ysq *(xnum + p(5)) / (xden + q(5)) - result = (sqrpi - result) / y - if (jintt /= 2) then - ysq = aint(y*sixten)/sixten - del = (y-ysq)*(y+ysq) - result = exp(-ysq*ysq) * exp(-del) * result - end if - else if (y >= xhuge) then - result = sqrpi / y - else - result = zero - end if -!------------------------------------------------------------------ -! fix up for negative argument, erf, etc. -!------------------------------------------------------------------ - if (jintt == 0) then - result = (half - result) + half - if (x < zero) result = -result - else if (jintt == 1) then - if (x < zero) result = two - result - else - if (x < zero) then - if (x < xneg) then - result = xinf - else - ysq = aint(x*sixten)/sixten - del = (x-ysq)*(x+ysq) - y = exp(ysq*ysq) * exp(del) - result = (y+y) - result - end if - end if - end if - end subroutine calerf - - function derf(x) -!-------------------------------------------------------------------- -! -! this subprogram computes approximate values for erf(x). -! (see comments heading calerf). -! -! author/date: w. j. cody, january 8, 1985 -! -!-------------------------------------------------------------------- - implicit none - real(wp_) :: derf - real(wp_), intent(in) :: x - integer :: jintt - real(wp_) :: result -!------------------------------------------------------------------ - jintt = 0 - call calerf(x,result,jintt) - derf = result - end function derf - - function derfc(x) -!-------------------------------------------------------------------- -! -! this subprogram computes approximate values for erfc(x). -! (see comments heading calerf). -! -! author/date: w. j. cody, january 8, 1985 -! -!-------------------------------------------------------------------- - implicit none - real(wp_) :: derfc - real(wp_), intent(in) :: x - integer :: jintt - real(wp_) :: result -!------------------------------------------------------------------ - jintt = 1 - call calerf(x,result,jintt) - derfc = result - end function derfc - - function derfcx(x) -!------------------------------------------------------------------ -! -! this subprogram computes approximate values for exp(x*x) * erfc(x). -! (see comments heading calerf). -! -! author/date: w. j. cody, march 30, 1987 -! -!------------------------------------------------------------------ - implicit none - real(wp_) :: derfcx - real(wp_), intent(in) :: x - integer :: jintt - real(wp_) :: result -!------------------------------------------------------------------ - jintt = 2 - call calerf(x,result,jintt) - derfcx = result - end function derfcx +! function derfcx(x) +!!------------------------------------------------------------------ +!! +!! this subprogram computes approximate values for exp(x*x) * erfc(x). +!! (see comments heading calerf). +!! +!! author/date: w. j. cody, march 30, 1987 +!! +!!------------------------------------------------------------------ +! implicit none +! real(wp_) :: derfcx +! real(wp_), intent(in) :: x +! integer :: jintt +! real(wp_) :: result +!!------------------------------------------------------------------ +! jintt = 2 +! call calerf(x,result,jintt) +! derfcx = result +! end function derfcx end module eierf \ No newline at end of file diff --git a/src/gray.f b/src/gray.f index 53abe33..1191d54 100644 --- a/src/gray.f +++ b/src/gray.f @@ -116,7 +116,7 @@ c c subroutine after_gray_integration use const_and_precisions, only : wp_,zero - use graydata_flags, only : ibeam,iwarm,ilarm,iequil,iprof, + use graydata_flags, only : ibeam,iwarm,iequil,iprof, . filenmeqq,filenmprf,filenmbm use graydata_anequil, only : dens0,te0 implicit none @@ -474,7 +474,7 @@ c c subroutine print_output(i,j,k) use const_and_precisions, only : wp_,pi - use graydata_flags, only : iequil,istpl0 + use graydata_flags, only : istpl0 use graydata_par, only : psdbnd implicit none c local constants @@ -673,7 +673,6 @@ c common/external functions/variables integer :: nstep,nrayr,nrayth real(wp_) :: xgcn,ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw, . phir,anx0c,any0c,anz0c,x00,y00,z00,bres,p0mw,sox,alpha0,beta0 - integer :: length c common/xgcn/xgcn common/nstep/nstep @@ -1116,7 +1115,6 @@ c local variables c common/external functions/variables real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,x00,y00,z00, . alpha0,beta0,ak0,akinv,fhz - integer :: length c common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 @@ -1221,10 +1219,10 @@ c use graydata_par, only : sgnbphi,sgniphi,factb use interp_eqprof, only : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz, . psia,psiant,psinop,btrcen,rcen,btaxis,rmaxis,zmaxis,rmnm, - . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,fpolas,phitedge,rrtor,rup,zup, + . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,fpolas,phitedge,rup,zup, . rlw,zlw,tr,tz,tfp,cc=>cceq,cfp,cc01=>cceq01,cc10=>cceq10, . cc20=>cceq20,cc02=>cceq02,cc11=>cceq11,psinr,qpsi,rv,zv,psin, - . psi,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd + . psi,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd!,rrtor implicit none c local constants integer, parameter :: kspl=3 @@ -2519,7 +2517,7 @@ c use const_and_precisions, only : wp_,pi use graydata_par, only : rwallm use magsurf_data, only : npoints - use interp_eqprof, only : psia,psiant,psinop,nsr=>nsrt,nsz=>nszt, + use interp_eqprof, only : psiant,psinop,nsr=>nsrt,nsz=>nszt, . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr use dierckx, only : profil,sproota implicit none @@ -3039,9 +3037,8 @@ c if(allocated(tjp)) deallocate(tjp) if(allocated(tlm)) deallocate(tlm) if(allocated(ch)) deallocate(ch) - if(allocated(ch01)) deallocate(ch01) allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), - . ch01(lw01),rctemp(npoints),zctemp(npoints),stat=ierr) + . rctemp(npoints),zctemp(npoints),stat=ierr) if (ierr.ne.0) return c c computation of flux surface averaged quantities @@ -3360,7 +3357,7 @@ c spline coefficients of rhot iopt=0 call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) -c spline interpolation of H(lambda,rhop) and dH/dlambda +c spline interpolation of H(lambda,rhop) iopt=0 s=0.0_wp_ @@ -3370,9 +3367,6 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda njpt=njp nlmt=nlm - call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, - . ch01,lw01,ier) -c return 99 format(20(1x,e12.5)) end @@ -5777,9 +5771,9 @@ c c subroutine pabs_curr(i,j,k) use const_and_precisions, only : wp_,pi - use graydata_flags, only : iequil,dst + use graydata_flags, only : dst use graydata_par, only : sgnbphi - use graydata_anequil, only : b0,rr0m,rpam + use graydata_anequil, only : rr0m implicit none c local constants integer, parameter :: jmx=31,kmx=36,nmx=8000 @@ -5937,20 +5931,24 @@ c . vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ use green_func_p, only: SpitzFuncCoeff use conical, only : fconic + use magsurf_data, only : ch,tjp,tlm,njpt,nlmt + use dierckx, only : profil implicit none c local constants real(wp_), parameter :: mc2m2=1.0_wp_/mc2**2, . canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi) + integer, parameter :: ksp=3 c arguments integer ieccd,nhmn,nhmx,ithn,iokhawa,ierr real(wp_) :: yg,anpl,anprre,amu,Zeff,rbn,rbx, . fc,dens,psinv,effjcd complex(wp_) :: ex,ey,ez c local variables - integer nhn + integer nhn,njp,nlm,npar real(wp_) :: anum,denom,resp,resj,coullog,anucc,alams, . fp0s,pa,cst2 - real(wp_), dimension(5) :: eccdpar + real(wp_) :: chlm(nlmt) + real(wp_), dimension(3+2*nlmt) :: eccdpar ! dimension(max(5,3+nlmt)) c real(wp_), external :: fjch,fjncl,fjch0 c @@ -6009,10 +6007,18 @@ c rzfc=(1+Zeff)/fc call SpitzFuncCoeff(amu,Zeff,fc) eccdpar(2) = fc eccdpar(3) = rbx - eccdpar(4) = psinv + + njp=njpt + nlm=nlmt + call profil(0,tjp,njp,tlm,nlm,ch,ksp,ksp,sqrt(psinv),nlm,chlm, + . ierr) + if(ierr.gt.0) print*,' Hlambda profil =',ierr + npar=3+2*nlm + eccdpar(4:3+nlm) = tlm + eccdpar(4+nlm:npar) = chlm do nhn=nhmn,nhmx call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, - . fjncl,eccdpar,4,resj,resp,iokhawa,ierr) + . fjncl,eccdpar,npar,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do @@ -6456,7 +6462,8 @@ c c extrapar(16) = zeff c extrapar(17) = fc c extrapar(18) = hb -c extrapar(19) = psin +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 @@ -6466,8 +6473,9 @@ c arguments real(wp_) :: upl,fjncl real(wp_), dimension(npar) :: extrapar c local variables - real(wp_) :: anpl,amu,ygn,zeff,fc,hb,psin,gam,u2,u,upr2, + real(wp_) :: anpl,amu,ygn,zeff,fc,hb,gam,u2,u,upr2, . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp + integer :: nlm c anpl=extrapar(2) amu=extrapar(3) @@ -6475,7 +6483,6 @@ c zeff=extrapar(16) fc=extrapar(17) hb=extrapar(18) - psin=extrapar(19) gam=anpl*upl+ygn u2=gam*gam-1.0_wp_ @@ -6483,12 +6490,15 @@ c upr2=u2-upl*upl bth=sqrt(2.0_wp_/amu) uth=u/bth - call GenSpitzFunc(amu,Zeff,fc,uth,u,gam,fk,dfk) + call GenSpitzFunc(Zeff,fc,uth,u,gam,fk,dfk) fk=fk*(4.0_wp_/amu**2) dfk=dfk*(2.0_wp_/amu)*bth alam=upr2/u2/hb - call vlambda(alam,psin,fu,dfu) + nlm=(npar-18)/2 + call vlambda(alam,extrapar(19), + . extrapar(19+nlm), + . nlm,fu,dfu) eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb if(upl.lt.0) eta=-eta @@ -6498,39 +6508,28 @@ c c c c - subroutine vlambda(alam,psi,fv,dfv) + subroutine vlambda(alam,tlm,chlm,nlmt,fv,dfv) use const_and_precisions, only : wp_ - use magsurf_data, only : ch,ch01,tjp,tlm,njpt,nlmt,npsi - use dierckx, only : fpbisp + use dierckx, only : splev,splder implicit none c local constants - integer, parameter :: nlam=41,ksp=3,nlest=nlam+ksp+1 + integer, parameter :: ksp=3 c arguments - real(wp_) alam,psi,fv,dfv + real(wp_) :: alam,fv,dfv + integer :: nlmt + real(wp_), dimension(nlmt) :: tlm,chlm c local variables - integer :: njp,nlm,iwp,iwl,njest,lwrk,kwrk - integer, dimension(:), allocatable :: iwrk - real(wp_), dimension(1) :: xxs,yys,ffs - real(wp_), dimension(:), allocatable :: wrk + integer :: nlm,ier + real(wp_), dimension(nlmt) :: wrk + real(wp_), dimension(1) :: xxs,ffs c - njest=npsi+ksp+1 - kwrk=npsi+nlam+njest+nlest+3 - lwrk=4*(npsi+nlam)+11*(njest+nlest)+njest*npsi+nlest+54 - allocate(iwrk(kwrk),wrk(lwrk)) -c - njp=njpt nlm=nlmt - xxs(1)=sqrt(psi) - yys(1)=alam + xxs(1)=alam c - call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs, - . wrk(1),wrk(5),iwrk(1),iwrk(2)) + call splev(tlm,nlm,chlm,ksp,xxs(1),ffs(1),1,ier) fv=ffs(1) c - iwp=1+(njp-4)*(nlm-5) - iwl=iwp+4 - call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2, - . xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2)) + call splder(tlm,nlm,chlm,ksp,1,xxs(1),ffs(1),1,wrk,ier) dfv=ffs(1) c return @@ -6631,8 +6630,7 @@ c subroutine pec(pabs,currt) use const_and_precisions, only : wp_,pi,one use numint, only : simpson - use graydata_flags, only : ipec,ieccd,iequil,nnd,dst - use graydata_anequil, only : rr0m,rpam + use graydata_flags, only : ipec,ieccd,nnd,dst use utils, only : locatex, locate, intlin implicit none c local constants @@ -6986,7 +6984,7 @@ c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) use const_and_precisions, only : wp_,emn1 - use graydata_flags, only : ipec,iequil + use graydata_flags, only : ipec use utils, only : locatex, locate, intlin, vmaxmini implicit none c arguments diff --git a/src/green_func_p.f90 b/src/green_func_p.f90 index 44eedbf..1a6f93a 100644 --- a/src/green_func_p.f90 +++ b/src/green_func_p.f90 @@ -13,7 +13,7 @@ IMPLICIT NONE CHARACTER(Len=1), PRIVATE :: adj_appr(6) ! adjoint approach switcher !------- - REAL(wp_), PRIVATE :: r2,q2,gp1,Rfactor + REAL(wp_), PRIVATE :: r2,q2,gp1 !------- REAL(wp_), PRIVATE, PARAMETER :: delta = 1e-4 ! border for recalculation !------- for N.M. subroutines (variational principle) ------- @@ -107,7 +107,7 @@ END SUBROUTINE Setup_SpitzFunc - SUBROUTINE GenSpitzFunc(mu,Zeff,fc,u,q,gam, K,dKdu) + SUBROUTINE GenSpitzFunc(Zeff,fc,u,q,gam, K,dKdu) !======================================================================= ! Author: N.B.Marushchenko ! June 2005: as start point the subroutine of Ugo Gasparino (198?) @@ -162,7 +162,6 @@ ! u - p/sqrt(2mT) ! q - p/mc; ! gam - relativistic factor; -! mu - mc2/Te ! Zeff - effective charge; ! fc - fraction of circulating particles. ! @@ -172,9 +171,9 @@ !======================================================================= use const_and_precisions, only : comp_eps IMPLICIT NONE - REAL(wp_), INTENT(in) :: mu,Zeff,fc,u,q,gam + REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam REAL(wp_), INTENT(out) :: K,dKdu - REAL(wp_) :: gam1,gam2,gam3,w,dwdu + REAL(wp_) :: gam1,gam2,gam3 !======================================================================= K = 0 dKdu = 0 @@ -241,7 +240,6 @@ INTEGER :: n,i,j REAL(wp_) :: rtc,rtc1,y,tn(1:nre) REAL(wp_) :: m(0:4,0:4),g(0:4) - REAL(wp_) :: om(0:4,0:4) REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, & gam22,gam32,gam42,gam02, & gam33,gam43,gam03, & @@ -250,9 +248,9 @@ alp23,alp24,alp20, & alp34,alp30,alp40 REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0 - LOGICAL :: renew,rel,newmu,newne,newZ,newfc + LOGICAL :: renew,rel,newmu,newZ,newfc REAL(wp_), SAVE :: sfdx(1:4) = 0 - REAL(wp_), SAVE :: ne_old =-1, mu_old =-1, Zeff_old =-1, fc_old =-1 + REAL(wp_), SAVE :: mu_old =-1, Zeff_old =-1, fc_old =-1 !======================================================================= rel = mu < mc2_ newmu = abs(mu -mu_old ) > delta*mu