simpson quadrature replaced with simpler rule for unevenly spaced intervals (computation of <rho>_{P,J} and d<rho>_{P,J)

src/gray-externals.f

@@ -6167,36 +6167,30 @@ c
         end do
       end do
  
-      h=1.0d0/dble(nd-1)    !!!!!!!!!! equi-spaced grid expected in
-      rhotpav=0.0d0         !!!!!!!!!! simpson subroutine
-      drhotpav=0.0d0        !!!!!!!!!! wrong integrals until fixed
-      rhotjav=0.0d0         !!!!!!!!!! (d)rhotpav,rhotjav,(d)rhotjava affected
+      rhotpav=0.0d0
+      rhot2pav=0.0d0
+      rhotjav=0.0d0
       rhotjava=0.0d0
       rhot2java=0.0d0
      
-      fi=dpdv/h
-      call simpson (nd,h,fi,spds)
-      fi=rhotv*dpdv/h
-      call simpson (nd,h,fi,rhotpav)
-      fi=rhotv*rhotv*dpdv/h
-      call simpson (nd,h,fi,rhot2pav)
-      rhotpav=rhotpav/spds
-      rhot2pav=rhot2pav/spds
+c     here dpdv is still dP (not divided yet by dV)
+      spds=sum(dpdv(1:nd))
+      if (spds.gt.0.0d0) then
+        rhotpav=sum(rhotv(1:nd)*dpdv(1:nd))/spds
+        rhot2pav=sum(rhotv(1:nd)*rhotv(1:nd)*dpdv(1:nd))/spds
+      end if
 
+c     here ajphiv is still dI (not divided yet by dA)
       if (ieccd.ne.0) then      
-      fi=ajphiv/h
-      call simpson (nd,h,fi,sccs)
-      fi=rhotv*ajphiv/h
-      call simpson (nd,h,fi,rhotjav)    
-      fi=abs(ajphiv)/h
-      call simpson (nd,h,fi,sccsa)
-      fi=rhotv*abs(ajphiv)/h
-      call simpson (nd,h,fi,rhotjava)
-      fi=rhotv*rhotv*abs(ajphiv)/h
-      call simpson (nd,h,fi,rhot2java)     
-      rhotjav=rhotjav/sccs
-      rhotjava=rhotjava/sccsa
-      rhot2java=rhot2java/sccsa
+        sccs=sum(ajphiv(1:nd))
+        if (sccs.gt.0.0d0) then
+          rhotjav=sum(rhotv(1:nd)*ajphiv(1:nd))/sccs
+        end if
+        sccsa=sum(abs(ajphiv(1:nd)))
+        if (sccsa.gt.0.0d0) then
+          rhotjava=sum(rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa
+          rhot2java=sum(rhotv(1:nd)*rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa
+        end if
       end if
 
 c     factor sqrt(8)=2 sqrt(2) to match with full width

src/grayl.f

@@ -257,6 +257,25 @@ c
       end subroutine simpson
 c
 c
+c
+      subroutine trapezoid(n,xi,fi,s)
+c     subroutine for integration with the trapezoidal rule.
+c     fi: integrand f(x); xi: abscissa x;
+c     s: integral Int_{xi(1)}^{xi(n)} f(x)dx
+        implicit none
+        integer n
+        real*8 xi(n),fi(n)
+        real*8 s
+        integer i
+c
+        s  = 0.0d0
+        do i = 1, n-1
+          s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
+        end do
+        s = 0.5d0*s
+      end subroutine trapezoid
+c
+c
 c     spline routines: begin
 c
 c