diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 4abc07a..a0d8744 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -145,7 +145,7 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) ! Yc = 1 -½ N∥² (weakly relativistic) ! if (weakly) then - Yc = max(one - npl2/2, zero) + Yc = max(1 - npl2/2, zero) else Yc = sqrt(max(1 - npl2, zero)) end if @@ -196,7 +196,7 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) else rdu2 = Yn**2 - Yc**2 gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2) - argexp = mu*(gamma - one) + argexp = mu*(gamma - 1) end if if (argexp <= expcr) then @@ -406,7 +406,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & Npr2 = Npr2_prev + 0.05_wp_ * (Npr2 - Npr2_prev) / abs(Npr2 - Npr2_prev) ! Again, make sure that we have a damped EM wave and not a ! Bernstein-like wave (see above) - if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < zero) & + if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < 0) & Npr2 = conjg(sqrt(Npr2))**2 end if end block modify_fixed_point @@ -464,7 +464,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & if (sox < 0) then e(3) = 1 else - e(1) = sqrt(one/(1 + abs(-eps(1,1)/eps(1,2))**2)) + e(1) = sqrt(1/(1 + abs(-eps(1,1)/eps(1,2))**2)) e(2) = -e(1)*eps(1,1)/eps(1,2) end if end if @@ -673,7 +673,7 @@ subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm) do n=-lrm,lrm do k=0,2 do m=0,lrm - rr(n,k,m)=zero + rr(n,k,m) = 0 end do end do end do @@ -893,7 +893,7 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm,error) do n=-lrm,lrm do k=0,2 do m=0,lrm - rr(n,k,m)=zero + rr(n,k,m) = 0 end do end do end do @@ -1059,7 +1059,7 @@ function fhermit(t,apar,npar) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) - ffe=zero + ffe = 0 uplh=upl**ih if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 if(m.eq.1) ffe=(one+s*(one-zm*fe0m))*uplh/mu2 @@ -1089,7 +1089,7 @@ subroutine antihermitian(ri,yg,mu,npl,ci,lrm) do n=1,lrm do k=0,2 do m=1,lrm - ri(n,k,m)=zero + ri(n,k,m) = 0 end do end do end do @@ -1101,7 +1101,7 @@ subroutine antihermitian(ri,yg,mu,npl,ci,lrm) do n=1,lrm ygn=n*yg rdu2=ygn**2-dnl - if(rdu2.gt.zero) then + if(rdu2.gt.0) then rdu=sqrt(rdu2) du=rdu/dnl ub=npl*ygn/dnl @@ -1231,17 +1231,17 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error) phim=sqrt(abs(phi2)) if (alpha.ge.0) then xp=psi-phim - yp=zero + yp=0 xm=-psi-phim - ym=zero + ym=0 x0=-phim - y0=zero + y0=0 else xp=psi yp=phim xm=-psi ym=phim - x0=zero + x0=0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) @@ -1295,19 +1295,19 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error) alpha=npl*npl/2.0_wp_-is*yg-one phi2=mu*alpha phim=sqrt(abs(phi2)) - if (alpha.ge.zero) then + if (alpha.ge.0) then xp=psi-phim - yp=zero + yp=0 xm=-psi-phim - ym=zero + ym=0 x0=-phim - y0=zero + y0=0 else xp=psi yp=phim xm=-psi ym=phim - x0=zero + x0=0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) @@ -1320,7 +1320,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error) ! cf12=czero if (alpha.ge.0) then - if (alpha.ne.zero) cf12=-(czp+czm)/(2.0_wp_*phim) + if (alpha.ne.0) cf12=-(czp+czm)/(2.0_wp_*phim) else cf12=-im*(czp+czm)/(2.0_wp_*phim) end if diff --git a/src/gray_core.f90 b/src/gray_core.f90 index bdfdde8..b0ac49e 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -131,10 +131,10 @@ contains allocate(results%jcd(params%output%nrho)) ! ...and initialise them - results%pabs = zero - results%icd = zero - results%dpdv = zero - results%jcd = zero + results%pabs = 0 + results%icd = 0 + results%dpdv = 0 + results%jcd = 0 ! ========= set environment END ========= ! Pre-determinted tables @@ -180,15 +180,15 @@ contains write (msg, '("pass: ",g0)') ip call log_info(msg, mod='gray_core', proc='gray_main') - pabs_pass = zero - icd_pass = zero + pabs_pass = 0 + icd_pass = 0 istop_pass = 0 ! stop flag for current pass nbeam_pass = 2*nbeam_pass ! max n of beams in current pass if(ip > 1) then - du1 = zero - gri = zero - ggri = zero + du1 = 0 + gri = 0 + ggri = 0 if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes end if @@ -215,16 +215,23 @@ contains cycle end if - call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) + psjki = 0 + ppabs = 0 + ccci = 0 + tau0 = 0 + alphaabs0 = 0 + dids0 = 0 + ccci0 = 0 + iiv = 1 if(ip == 1) then ! 1st pass igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass - tau1 = zero ! * tau from previous passes - etau1 = one - cpl1 = one ! * coupling from previous passes - lgcpl1 = zero - p0ray = p0jk ! * initial beam power + tau1 = 0 ! * tau from previous passes + etau1 = 1 + cpl1 = 1 ! * coupling from previous passes + lgcpl1 = 0 + p0ray = p0jk ! * initial beam power call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri) @@ -469,10 +476,10 @@ contains end if end block else - tekev=zero - alpha=zero - didp=zero - anprim=zero + tekev=0 + alpha=0 + didp=0 + anprim=0 anprre=anpr nharm=0 nhf=0 @@ -568,15 +575,15 @@ contains cpl_beam1 = sum(& p0ray * exp(-tau0) * cpls(:,child_index_rt)/cpl1, MASK=iop > 2) / & sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) - cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam + cpl_beam2 = 1 - cpl_beam1 ! * average X-mode coupling for next beam if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1) - cpl_cbeam2 = one - cpl_cbeam1 + cpl_cbeam2 = 1 - cpl_cbeam1 end if else ! last pass OR no ray re-entered plasma - cpl_beam1 = zero - cpl_beam2 = zero + cpl_beam1 = 0 + cpl_beam2 = 0 end if ! print final results for pass on screen @@ -661,36 +668,6 @@ contains end subroutine gray_main - - - subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) - use const_and_precisions, only : zero -! arguments - real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci - real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0 - integer, dimension(:), intent(out) :: iiv -!! common/external functions/variables -! integer :: jclosest -! real(wp_), dimension(3) :: anwcl,xwcl -! -! common/refln/anwcl,xwcl,jclosest -! -! jclosest=nrayr+1 -! anwcl(1:3)=0.0_wp_ -! xwcl(1:3)=0.0_wp_ - - psjki = zero - ppabs = zero - ccci = zero - tau0 = zero - alphaabs0 = zero - dids0 = zero - ccci0 = zero - iiv = 1 - - end subroutine vectinit - - subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, & pos, grad_u, grad, hess) ! Computes the initial conditions for tracing a beam @@ -1007,7 +984,6 @@ contains subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error) - use const_and_precisions, only : zero, half use gray_params, only : raytracing_parameters use gray_errors, only : gray_error, unstable_beam, raise_error @@ -1060,7 +1036,7 @@ contains call solg0(dxv1,dxv2,dxv3,dgu) du1(:,k,j) = dgu end do - gri(:,1) = zero + gri(:,1) = 0 ! compute grad u1 and grad(S_I) for all the other rays if (params%nrayr > 1) then @@ -1125,7 +1101,7 @@ contains end do ! compute derivatives of grad u and grad(S_I) for rays jk>1 - ggri(:,:,1) = zero + ggri(:,:,1) = 0 jm=1 j=2 k=0 @@ -1160,9 +1136,9 @@ contains uxx = dgg(1,1) uyy = dgg(2,2) uzz = dgg(3,3) - uxy = (dgg(1,2) + dgg(2,1))*half - uxz = (dgg(1,3) + dgg(3,1))*half - uyz = (dgg(2,3) + dgg(3,2))*half + uxy = (dgg(1,2) + dgg(2,1))/2 + uxz = (dgg(1,3) + dgg(3,1))/2 + uyz = (dgg(2,3) + dgg(3,2))/2 ! derivatives of S_I and Grad(S_I) gx = ux*dffiu @@ -1243,7 +1219,7 @@ contains subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg, psinv_opt) - use const_and_precisions, only : zero, cm + use const_and_precisions, only : cm use gray_equil, only : abstract_equil, vacuum use gray_plasma, only : abstract_plasma @@ -1267,15 +1243,15 @@ contains real(wp_) :: brr,bphi,bzz,dxgdpsi real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi - xg = zero + xg = 0 yg = 99._wp_ psinv = -1._wp_ - dens = zero - btot = zero - derxg = zero - deryg = zero - bv = zero - derbv = zero + dens = 0 + btot = 0 + derxg = 0 + deryg = 0 + bv = 0 + derbv = 0 select type (equil) type is (vacuum) @@ -1284,11 +1260,11 @@ contains return end select - dbtot = zero - dbv = zero - dbvcdc = zero - dbvcdc = zero - dbvdc = zero + dbtot = 0 + dbv = 0 + dbvcdc = 0 + dbvcdc = 0 + dbvdc = 0 xx = xv(1) yy = xv(2) @@ -1315,7 +1291,7 @@ contains if (present(psinv_opt)) psinv_opt = psinv ! compute yg and derivative - if(psinv < zero) then + if(psinv < 0) then bphi = fpolv/rrm btot = abs(bphi) yg = btot/bres @@ -1405,7 +1381,7 @@ contains ! (`ywppla_upd` suborutine); while the optional ones are used for ! computing the absoprtion and current drive. - use const_and_precisions, only : zero, one, half, two + use const_and_precisions, only : zero, half use gray_params, only : gray_parameters, STEP_ARCLEN, & STEP_TIME, STEP_PHASE @@ -1459,26 +1435,26 @@ contains anpl = dot_product(anv, bv) ! N∥ = N̅⋅B̅ ! Shorthands used in the expressions below - yg2 = yg**2 ! Y² - anpl2 = anpl**2 ! N∥² - dnl = one - anpl2 ! 1 - N∥² - duh = one - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance) + yg2 = yg**2 ! Y² + anpl2 = anpl**2 ! N∥² + dnl = 1 - anpl2 ! 1 - N∥² + duh = 1 - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance) ! Compute/copy optional outputs if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥ if (present(anpl_)) anpl_ = anpl ! N∥ - an2s = one - dan2sdxg = zero - dan2sdyg = zero - dan2sdnpl = zero - del = zero - fdia = zero - dfdiadnpl = zero - dfdiadxg = zero - dfdiadyg = zero + an2s = 1 + dan2sdxg = 0 + dan2sdyg = 0 + dan2sdnpl = 0 + del = 0 + fdia = 0 + dfdiadnpl = 0 + dfdiadxg = 0 + dfdiadyg = 0 - if(xg > zero) then + if(xg > 0) then ! Derivatives of the cold plasma refractive index ! ! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)] @@ -1487,19 +1463,19 @@ contains ! + for the X mode, - for the O mode ! √Δ - del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) + del = sqrt(dnl**2 + 4.0_wp_*anpl2*(1 - xg)/yg2) ! ∂(N²s)/∂X ! Note: this term is nonzero for X=0, but it multiplies terms ! proportional to X or ∂X/∂ψ which are zero outside the plasma. - dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 & - + sox*xg*anpl2/(del*duh) - one + dan2sdxg = - half*yg2*(1 - yg2)*(1 + anpl2 + sox*del)/duh**2 & + + sox*xg*anpl2/(del*duh) - 1 ! ∂(N²s)/∂Y - dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 & - + two*sox*xg*(one - xg)*anpl2/(yg*del*duh) + dan2sdyg = - xg*yg*(1 - xg)*(1 + anpl2 + sox*del)/duh**2 & + + 2*sox*xg*(1 - xg)*anpl2/(yg*del*duh) ! ∂(N²s)/∂N∥ - dan2sdnpl = - xg*yg2*anpl/duh & - - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) + dan2sdnpl = - xg*yg2*anpl/duh & + - sox*xg*anpl*(2*(1 - xg) - yg2*dnl)/(del*duh) if(igrad > 0) then ! Derivatives used in the complex eikonal terms (beamtracing only) @@ -1507,28 +1483,27 @@ contains real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel ! ∂²Δ/∂N∥² - ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & + ddelnpl2 = 2*(2*(1 - xg)*(1 + 3.0_wp_*anpl2**2) & - yg2*dnl**3)/yg2/del**3 ! ∂²(N²s)/∂N∥² - fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh + fdia = - xg*yg2*(1 + half*sox*ddelnpl2)/duh ! Intermediates results used right below - derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & - - dnl**2*(one + 3.0_wp_*anpl2)*yg2 - derdel = 4.0_wp_*derdel/(yg*del)**5 - ddelnpl2y = two*(one - xg)*derdel + derdel = 2*(1 - xg)*anpl2*(1 + 3*anpl2**2) & + - dnl**2*(1 + 3*anpl2)*yg2 + derdel = 4*derdel/(yg*del)**5 + ddelnpl2y = 2*(1 - xg)*derdel ddelnpl2x = yg*derdel ! ∂³(N²s)/∂N∥³ - dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & - /(yg2*del**5) + dfdiadnpl = 24*sox*xg*(1 - xg)*anpl*(1 - anpl2**2)/(yg2*del**5) ! ∂³(N²s)/∂N∥²∂X - dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & - *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) + dfdiadxg = - yg2*(1 - yg2)/duh**2 - sox*yg2*((1 - yg2) & + *ddelnpl2 + xg*duh*ddelnpl2x)/(2*duh**2) ! ∂³(N²s)/∂N∥²∂Y - dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & - - sox*xg*yg*(two*(one - xg)*ddelnpl2 & - + yg*duh*ddelnpl2y)/(two*duh**2) + dfdiadyg = - 2*yg*xg*(1 - xg)/duh**2 & + - sox*xg*yg*(2*(1 - xg)*ddelnpl2 & + + yg*duh*ddelnpl2y)/(2*duh**2) end block end if end if @@ -1541,7 +1516,7 @@ contains ! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅ ! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅. - derdnv = two*anv - dan2sdnpl*bv + derdnv = 2*anv - dan2sdnpl*bv ! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω ! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω @@ -1549,7 +1524,7 @@ contains ! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason. ! 3. N gains a dependency on ω because Λ(∇S, ω) is computed ! on the constrains Λ=0. - derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl + derdom = -2*an2 + 2*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl if (igrad > 0) then ! Complex eikonal terms added to the above expressions @@ -1576,7 +1551,7 @@ contains ! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²) ! Note: as above ∇S_I gains a dependency on ω - derdom = derdom + two*gr2 - bdotgr**2 & + derdom = derdom + 2*gr2 - bdotgr**2 & * (fdia + xg*dfdiadxg + half*yg*dfdiadyg & + half*anpl*dfdiadnpl) end block @@ -1651,7 +1626,7 @@ contains if (present(ddr) .or. present(ddi)) then ! Dispersion relation (geometric optics) ! ddr ← Λ = N² - N²s(X,Y,N∥) = 0 - an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh + an2s = 1 - xg - half*xg*yg2*(1 + anpl2 + sox*del)/duh ddr = an2 - an2s end if diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90 index e5e799f..30a198e 100644 --- a/src/gray_plasma.f90 +++ b/src/gray_plasma.f90 @@ -185,7 +185,7 @@ contains ddens = self%dens_spline%deriv(psin) ! Evaluate the spline 1st derivative - if (abs(dens) < 1.0e-10_wp_) dens = zero + if (abs(dens) < 1.0e-10_wp_) dens = 0 else ! Use a C² polynomial extension outside (ψ > ψ₀) diff --git a/src/reflections.f90 b/src/reflections.f90 index 153f896..9459266 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -1,5 +1,5 @@ module reflections - use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one + use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge implicit none @@ -16,7 +16,7 @@ subroutine reflect(ki,nsurf,ko) real(wp_), intent(out), dimension(3) :: ko real(wp_) :: twokn,norm2 norm2 = dot_product(nsurf,nsurf) - if (norm2>zero) then + if (norm2>0) then twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2 ko=ki-twokn*nsurf else @@ -42,20 +42,20 @@ subroutine inters_linewall(xv, kv, wall, sint, normw) sint=comp_huge iint=0 - normw=zero + normw=0 do i=1, size(wall%R)-1 ! search intersections with i-th wall segment call linecone_coord(xv, kv, wall%R(i:i+1), wall%z(i:i+1), si, ti, ni) ! discard solutions with s<=0 first=ni+1 do j=1,ni - if (si(j)>zero) then + if (si(j)>0) then first=j exit end if end do do j=first,ni - if ((si(j)=zero .and. ti(j)<=one) then + if ((si(j)=0 .and. ti(j)<=1) then ! check intersection is in r,z range and keep the closest sint = si(j) iint = i @@ -72,7 +72,7 @@ subroutine inters_linewall(xv, kv, wall, sint, normw) l = hypot(drw, dzw) kxy = norm2(kv(1:2)) normw(3) = -drw/l - if (rint>zero) then + if (rint>0) then normw(1) = xint/rint*dzw/l normw(2) = yint/rint*dzw/l else @@ -80,7 +80,7 @@ subroutine inters_linewall(xv, kv, wall, sint, normw) normw(2) = kv(2)/kxy*dzw/l end if ! reverse normal if k.n>0 - if (dot_product(normw,kv)>zero) normw=-normw + if (dot_product(normw,kv)>0) normw=-normw end subroutine inters_linewall @@ -168,8 +168,8 @@ subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr) dyb = yb(2)-yb(1) crossprod = dxb*dya - dxa*dyb if (abs(crossprod)=zero .and. s<=one .and. & - t>=zero .and. t<=one) interssegm = .true. + if (ierr==0 .and. s>=0 .and. s<=1 .and. & + t>=0 .and. t<=1) interssegm = .true. end function interssegm diff --git a/src/vendor/eierf.f90 b/src/vendor/eierf.f90 index 716f629..a5d481c 100644 --- a/src/vendor/eierf.f90 +++ b/src/vendor/eierf.f90 @@ -180,7 +180,7 @@ contains 1.99999999999048104167_wp_/) !---------------------------------------------------------------------- x = arg - if (x == zero) then + if (x == 0) then ei = -xinf if (intt == 2) ei = -ei else if ((x < zero) .or. (intt == 2)) then @@ -209,7 +209,7 @@ contains if (intt /= 3) ei = ei * exp(-y) else if ((y > xbig) .and. (intt < 3)) then - ei = zero + ei = 0 else w = one / y sump = e(1) @@ -231,8 +231,8 @@ contains !---------------------------------------------------------------------- t = x + x t = t / three - two - px(1) = zero - qx(1) = zero + px(1) = 0 + qx(1) = 0 px(2) = p(1) qx(2) = q(1) do i = 2, 9 @@ -262,14 +262,14 @@ contains if (intt == 3) ei = exp(-x) * ei end if else if (x < twelve) then - frac = zero + frac = 0 do i = 1, 9 frac = s(i) / (r(i) + x + frac) end do ei = (r(10) + frac) / x if (intt /= 3) ei = ei * exp(x) else if (x <= two4) then - frac = zero + frac = 0 do i = 1, 9 frac = q1(i) / (p1(i) + x + frac) end do @@ -280,7 +280,7 @@ contains ei = xinf else y = one / x - frac = zero + frac = 0 do i = 1, 9 frac = q2(i) / (p2(i) + x + frac) end do @@ -532,7 +532,7 @@ contains 1.99999999999048104167_wp_/) !---------------------------------------------------------------------- x = arg - if (x == zero) then + if (x == 0) then ei = -xinf else if ((x < zero)) then !---------------------------------------------------------------------- @@ -574,8 +574,8 @@ contains !---------------------------------------------------------------------- t = x + x t = t / three - two - px(1) = zero - qx(1) = zero + px(1) = 0 + qx(1) = 0 px(2) = p(1) qx(2) = q(1) do i = 2, 9 @@ -603,20 +603,20 @@ contains ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0 end if else if (x < twelve) then - frac = zero + frac = 0 do i = 1, 9 frac = s(i) / (r(i) + x + frac) end do ei = (r(10) + frac) / x else if (x <= two4) then - frac = zero + frac = 0 do i = 1, 9 frac = q1(i) / (p1(i) + x + frac) end do ei = (p1(10) + frac) / x else y = one / x - frac = zero + frac = 0 do i = 1, 9 frac = q2(i) / (p2(i) + x + frac) end do diff --git a/src/vendor/minpack.f90 b/src/vendor/minpack.f90 index e5a65d8..9217aa4 100644 --- a/src/vendor/minpack.f90 +++ b/src/vendor/minpack.f90 @@ -48,7 +48,6 @@ module minpack contains pure subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa) - use const_and_precisions, only : zero, one ! arguments integer, intent(in) :: n, ldfjac, lwa integer, intent(out) :: info @@ -166,7 +165,7 @@ contains ! ! check the input parameters for errors. ! - if (n <= 0 .or. ldfjac < n .or. tol < zero & + if (n <= 0 .or. ldfjac < n .or. tol < 0 & .or. lwa < (n*(n + 13))/2) return ! ! call hybrj. @@ -175,7 +174,7 @@ contains xtol = tol mode = 2 do j = 1, n - wa(j) = one + wa(j) = 1 end do nprint = 0 lr = (n*(n + 1))/2 @@ -188,7 +187,7 @@ contains pure subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, & factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, & wa3,wa4) - use const_and_precisions, only : zero, one, epsmch=>comp_eps + use const_and_precisions, only : epsmch=>comp_eps use, intrinsic :: ieee_exceptions, only : ieee_get_halting_mode, & ieee_set_halting_mode, & ieee_invalid @@ -373,12 +372,12 @@ contains ! ! check the input parameters for errors. ! - if (n <= 0 .or. ldfjac < n .or. xtol < zero & - .or. maxfev <= 0 .or. factor <= zero & + if (n <= 0 .or. ldfjac < n .or. xtol < 0 & + .or. maxfev <= 0 .or. factor <= 0 & .or. lr < (n*(n + 1))/2) go to 300 if (mode == 2) then do j = 1, n - if (diag(j) <= zero) go to 300 + if (diag(j) <= 0) go to 300 end do end if ! @@ -422,7 +421,7 @@ contains if (mode /= 2) then do j = 1, n diag(j) = wa2(j) - if (wa2(j) == zero) diag(j) = one + if (wa2(j) == 0) diag(j) = 1 end do end if ! @@ -434,7 +433,7 @@ contains end do xnorm = enorm(n,wa3) delta = factor*xnorm - if (delta == zero) delta = factor + if (delta == 0) delta = factor end if ! ! form (q transpose)*fvec and store in qtf. @@ -443,8 +442,8 @@ contains qtf(i) = fvec(i) end do do j = 1, n - if (fjac(j,j) /= zero) then - summ = zero + if (fjac(j,j) /= 0) then + summ = 0 do i = j, n summ = summ + fjac(i,j)*qtf(i) end do @@ -466,7 +465,7 @@ contains l = l + n - i end do r(l) = wa1(j) - if (wa1(j) == zero) sing = .true. + if (wa1(j) == 0) sing = .true. end do ! ! accumulate the orthogonal factor in fjac. @@ -520,14 +519,14 @@ contains ! ! compute the scaled actual reduction. ! - actred = -one - if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2 + actred = -1 + if (fnorm1 < fnorm) actred = 1 - (fnorm1/fnorm)**2 ! ! compute the scaled predicted reduction. ! l = 1 do i = 1, n - summ = zero + summ = 0 do j = i, n summ = summ + r(l)*wa1(j) l = l + 1 @@ -535,14 +534,14 @@ contains wa3(i) = qtf(i) + summ end do temp = enorm(n,wa3) - prered = zero - if (temp < fnorm) prered = one - (temp/fnorm)**2 + prered = 0 + if (temp < fnorm) prered = 1 - (temp/fnorm)**2 ! ! compute the ratio of the actual to the predicted ! reduction. ! - ratio = zero - if (prered > zero) ratio = actred/prered + ratio = 0 + if (prered > 0) ratio = actred/prered ! ! update the step bound. ! @@ -554,7 +553,7 @@ contains ncfail = 0 ncsuc = ncsuc + 1 if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5) - if (abs(ratio-one) <= p1) delta = pnorm/p5 + if (abs(ratio-1) <= p1) delta = pnorm/p5 end if ! ! test for successful iteration. @@ -582,7 +581,7 @@ contains ! ! test for convergence. ! - if (delta <= xtol*xnorm .or. fnorm == zero) info = 1 + if (delta <= xtol*xnorm .or. fnorm == 0) info = 1 if (info /= 0) go to 300 ! ! tests for termination and stringent tolerances. @@ -601,7 +600,7 @@ contains ! and update qtf if necessary. ! do j = 1, n - summ = zero + summ = 0 do i = 1, n summ = summ + fjac(i,j)*wa4(i) end do @@ -638,7 +637,7 @@ contains end subroutine hybrj pure subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2) - use const_and_precisions, only : zero, one, epsmch=>comp_eps + use const_and_precisions, only : epsmch=>comp_eps ! arguments integer, intent(in) :: n, lr real(wp_), intent(in) :: delta, r(lr), diag(n), qtb(n) @@ -712,20 +711,20 @@ contains jp1 = j + 1 jj = jj - k l = jj + 1 - summ = zero + summ = 0 do i = jp1, n summ = summ + r(l)*x(i) l = l + 1 end do temp = r(jj) - if (temp == zero) then + if (temp == 0) then l = j do i = 1, j temp = dmax1(temp,abs(r(l))) l = l + n - i end do temp = epsmch*temp - if (temp == zero) temp = epsmch + if (temp == 0) temp = epsmch end if x(j) = (qtb(j) - summ)/temp end do @@ -733,7 +732,7 @@ contains ! test whether the gauss-newton direction is acceptable. ! do j = 1, n - wa1(j) = zero + wa1(j) = 0 wa2(j) = diag(j)*x(j) end do qnorm = enorm(n,wa2) @@ -756,9 +755,9 @@ contains ! the special case in which the scaled gradient is zero. ! gnorm = enorm(n,wa1) - sgnorm = zero + sgnorm = 0 alpha = delta/qnorm - if (gnorm /= zero) then + if (gnorm /= 0) then ! ! calculate the point along the scaled gradient ! at which the quadratic is minimized. @@ -768,7 +767,7 @@ contains end do l = 1 do j = 1, n - summ = zero + summ = 0 do i = j, n summ = summ + r(l)*wa1(i) l = l + 1 @@ -780,7 +779,7 @@ contains ! ! test whether the scaled gradient direction is acceptable. ! - alpha = zero + alpha = 0 if (sgnorm < delta) then ! ! the scaled gradient direction is not acceptable. @@ -791,22 +790,21 @@ contains temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta) temp = temp - (delta/qnorm)*(sgnorm/delta)**2 & + sqrt((temp-(delta/qnorm))**2 & - +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2)) - alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp + +(1-(delta/qnorm)**2)*(1-(sgnorm/delta)**2)) + alpha = ((delta/qnorm)*(1 - (sgnorm/delta)**2))/temp end if end if ! ! form appropriate convex combination of the gauss-newton ! direction and the scaled gradient direction. ! - temp = (one - alpha)*dmin1(sgnorm,delta) + temp = (1 - alpha)*dmin1(sgnorm,delta) do j = 1, n x(j) = temp*wa1(j) + alpha*x(j) end do end subroutine dogleg pure function enorm(n,x) - use const_and_precisions, only : zero, one real(wp_) :: enorm integer, intent(in) :: n real(wp_), dimension(n), intent(in) :: x @@ -850,11 +848,11 @@ contains integer :: i real(wp_) :: agiant,floatn,s1,s2,s3,xabs,x1max,x3max real(wp_), parameter :: rdwarf=3.834e-20_wp_,rgiant=1.304e19_wp_ - s1 = zero - s2 = zero - s3 = zero - x1max = zero - x3max = zero + s1 = 0 + s2 = 0 + s3 = 0 + x1max = 0 + x3max = 0 floatn = n agiant = rgiant/floatn do i = 1, n @@ -865,7 +863,7 @@ contains ! sum for large components. ! if (xabs > x1max) then - s1 = one + s1*(x1max/xabs)**2 + s1 = 1 + s1*(x1max/xabs)**2 x1max = xabs else s1 = s1 + (xabs/x1max)**2 @@ -875,10 +873,10 @@ contains ! sum for small components. ! if (xabs > x3max) then - s3 = one + s3*(x3max/xabs)**2 + s3 = 1 + s3*(x3max/xabs)**2 x3max = xabs else - if (xabs /= zero) s3 = s3 + (xabs/x3max)**2 + if (xabs /= 0) s3 = s3 + (xabs/x3max)**2 end if end if else @@ -891,11 +889,11 @@ contains ! ! calculation of norm. ! - if (s1 /= zero) then + if (s1 /= 0) then enorm = x1max*sqrt(s1+(s2/x1max)/x1max) else - if (s2 /= zero) then - if (s2 >= x3max) enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3))) + if (s2 /= 0) then + if (s2 >= x3max) enorm = sqrt(s2*(1+(x3max/s2)*(x3max*s3))) if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max)+(x3max*s3))) else enorm = x3max*sqrt(s3) @@ -904,7 +902,6 @@ contains end function enorm pure subroutine qform(m,n,q,ldq,wa) - use const_and_precisions, only : zero, one ! arguments integer, intent(in) :: m,n,ldq real(wp_), intent(out) :: wa(m) @@ -956,7 +953,7 @@ contains do j = 2, minmn jm1 = j - 1 do i = 1, jm1 - q(i,j) = zero + q(i,j) = 0 end do end do ! @@ -965,9 +962,9 @@ contains np1 = n + 1 do j = np1, m do i = 1, m - q(i,j) = zero + q(i,j) = 0 end do - q(j,j) = one + q(j,j) = 1 end do ! ! accumulate q from its factored form. @@ -976,12 +973,12 @@ contains k = minmn - l + 1 do i = k, m wa(i) = q(i,k) - q(i,k) = zero + q(i,k) = 0 end do - q(k,k) = one - if (wa(k) /= zero) then + q(k,k) = 1 + if (wa(k) /= 0) then do j = k, m - summ = zero + summ = 0 do i = k, m summ = summ + q(i,j)*wa(i) end do @@ -995,7 +992,7 @@ contains end subroutine qform pure subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) - use const_and_precisions, only : zero, one, epsmch=>comp_eps + use const_and_precisions, only : zero, epsmch=>comp_eps ! arguments integer, intent(in) :: m, n, lda, lipvt integer, intent(out) :: ipvt(lipvt) @@ -1122,19 +1119,19 @@ contains ! j-th column of a to a multiple of the j-th unit vector. ! ajnorm = enorm(m-j+1,a(j,j)) - if (ajnorm /= zero) then - if (a(j,j) < zero) ajnorm = -ajnorm + if (ajnorm /= 0) then + if (a(j,j) < 0) ajnorm = -ajnorm do i = j, m a(i,j) = a(i,j)/ajnorm end do - a(j,j) = a(j,j) + one + a(j,j) = a(j,j) + 1 ! ! apply the transformation to the remaining columns ! and update the norms. ! jp1 = j + 1 do k = jp1, n - summ = zero + summ = 0 do i = j, m summ = summ + a(i,j)*a(i,k) end do @@ -1142,9 +1139,9 @@ contains do i = j, m a(i,k) = a(i,k) - temp*a(i,j) end do - if (pivot .and. rdiag(k) /= zero) then + if (pivot .and. rdiag(k) /= 0) then temp = a(j,k)/rdiag(k) - rdiag(k) = rdiag(k)*sqrt(dmax1(zero,one-temp**2)) + rdiag(k) = rdiag(k)*sqrt(dmax1(zero, 1-temp**2)) if (p05*(rdiag(k)/wa(k))**2 <= epsmch) then rdiag(k) = enorm(m-j,a(jp1,k)) wa(k) = rdiag(k) @@ -1157,7 +1154,6 @@ contains end subroutine qrfac pure subroutine r1mpyq(m,n,a,lda,v,w) - use const_and_precisions, only : one ! arguments integer, intent(in) :: m, n, lda real(wp_), intent(in) :: v(n),w(n) @@ -1221,10 +1217,10 @@ contains if (nm1 < 1) return do nmj = 1, nm1 j = n - nmj - if (abs(v(j)) > one) cs = one/v(j) - if (abs(v(j)) > one) sn = sqrt(one-cs**2) - if (abs(v(j)) <= one) sn = v(j) - if (abs(v(j)) <= one) cs = sqrt(one-sn**2) + if (abs(v(j)) > 1) cs = 1/v(j) + if (abs(v(j)) > 1) sn = sqrt(1-cs**2) + if (abs(v(j)) <= 1) sn = v(j) + if (abs(v(j)) <= 1) cs = sqrt(1-sn**2) do i = 1, m temp = cs*a(i,j) - sn*a(i,n) a(i,n) = sn*a(i,j) + cs*a(i,n) @@ -1235,10 +1231,10 @@ contains ! apply the second set of givens rotations to a. ! do j = 1, nm1 - if (abs(w(j)) > one) cs = one/w(j) - if (abs(w(j)) > one) sn = sqrt(one-cs**2) - if (abs(w(j)) <= one) sn = w(j) - if (abs(w(j)) <= one) cs = sqrt(one-sn**2) + if (abs(w(j)) > 1) cs = 1/w(j) + if (abs(w(j)) > 1) sn = sqrt(1-cs**2) + if (abs(w(j)) <= 1) sn = w(j) + if (abs(w(j)) <= 1) cs = sqrt(1-sn**2) do i = 1, m temp = cs*a(i,j) + sn*a(i,n) a(i,n) = -sn*a(i,j) + cs*a(i,n) @@ -1248,7 +1244,7 @@ contains end subroutine r1mpyq pure subroutine r1updt(m,n,s,ls,u,v,w,sing) - use const_and_precisions, only : zero, one, giant=>comp_huge + use const_and_precisions, only : giant=>comp_huge ! arguments integer, intent(in) :: m, n, ls logical, intent(out) :: sing @@ -1346,8 +1342,8 @@ contains do nmj = 1, nm1 j = n - nmj jj = jj - (m - j + 1) - w(j) = zero - if (v(j) /= zero) then + w(j) = 0 + if (v(j) /= 0) then ! ! determine a givens rotation which eliminates the ! j-th element of v. @@ -1356,8 +1352,8 @@ contains cotan = v(n)/v(j) sn = p5/sqrt(p25+p25*cotan**2) cs = sn*cotan - tau = one - if (abs(cs)*giant > one) tau = one/cs + tau = 1 + if (abs(cs)*giant > 1) tau = 1/cs else tn = v(j)/v(n) cs = p5/sqrt(p25+p25*tn**2) @@ -1393,7 +1389,7 @@ contains ! sing = .false. do j = 1, nm1 - if (w(j) /= zero) then + if (w(j) /= 0) then ! ! determine a givens rotation which eliminates the ! j-th element of the spike. @@ -1402,8 +1398,8 @@ contains cotan = s(jj)/w(j) sn = p5/sqrt(p25+p25*cotan**2) cs = sn*cotan - tau = one - if (abs(cs)*giant > one) tau = one/cs + tau = 1 + if (abs(cs)*giant > 1) tau = 1/cs else tn = w(j)/s(jj) cs = p5/sqrt(p25+p25*tn**2) @@ -1429,7 +1425,7 @@ contains ! ! test for zero diagonal elements in the output s. ! - if (s(jj) == zero) sing = .true. + if (s(jj) == 0) sing = .true. jj = jj + (m - j + 1) end do ! @@ -1440,7 +1436,7 @@ contains s(l) = w(i) l = l + 1 end do - if (s(jj) == zero) sing = .true. + if (s(jj) == 0) sing = .true. ! end subroutine r1updt diff --git a/src/vendor/numint.f90 b/src/vendor/numint.f90 index 036ac2a..3bdb63c 100644 --- a/src/vendor/numint.f90 +++ b/src/vendor/numint.f90 @@ -1,6 +1,6 @@ module numint - use const_and_precisions, only : wp_, zero, one + use const_and_precisions, only : wp_ implicit none @@ -16,10 +16,10 @@ contains integer :: i real(wp_) :: s0,s1,s2 - s = zero - s0 = zero - s1 = zero - s2 = zero + s = 0 + s0 = 0 + s1 = 0 + s2 = 0 do i = 2, n-1, 2 s1 = s1+fi(i-1) s0 = s0+fi(i) @@ -39,7 +39,7 @@ contains real(wp_), intent(out) :: s integer :: i - s = zero + s = 0 do i = 1, n-1 s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i)) end do @@ -110,11 +110,11 @@ contains ! ! initialize running sums to zero. ! - flag = zero - result = zero - cor11 = zero - errest = zero - area = zero + flag = 0 + result = 0 + cor11 = 0 + errest = 0 + area = 0 nofun = 0 if (a .eq. b) return ! @@ -124,7 +124,7 @@ contains nim = 1 x0 = a x(16) = b - qprev = zero + qprev = 0 f0 = fun(x0) stone = (b - a) / 16.0_wp_ x(8) = (x0 + x(16)) / 2.0_wp_ @@ -174,7 +174,7 @@ contains ! ! current level is levmax. ! - flag = flag + one + flag = flag + 1 exit end if if (nofun .gt. nofin) then @@ -243,7 +243,7 @@ contains ! ! make sure errest not less than roundoff level. ! - if (errest .eq. zero) return + if (errest .eq. 0) return do temp = abs(result) + errest if (temp .ne. abs(result)) return