Module coquad1
  Use quad1
  Implicit None
  Private co_rectangle_dp
  Public co_rectangle_rule
  Interface co_rectangle_rule
    Module Procedure co_rectangle_dp
  End Interface
Contains
  Subroutine co_rectangle_dp(f,dfrom,dto,dnstep,integral)
    Interface
      Double Precision Function f(x)
        Double Precision,Intent(In) :: x
      End Function
    End Interface
    Double Precision,Intent(In) :: dfrom,dto
    Integer,Intent(In) :: dnstep
    Double Precision,Intent(Out) :: integral[*]
    Double Precision from,to
    Integer i,nstep
    ! Calculate steps and bounds for this image.
    nstep = dnstep/Num_Images()
    from = (dto*(This_Image()-1)+dfrom*(Num_Images()-This_Image()+1))/Num_Images()
    to = (dto*(This_Image())+dfrom*(Num_Images()-This_Image()))/Num_Images()
    ! Use the single-image code to calculate the integral on our part of the interval.
    Call rectangle_rule(f,from,to,nstep,integral)
    ! We have integrated our interval, now synchronise.
    Sync All
    If (This_Image()==1) Then
      ! Everyone has finished, so sum all the values.
      Do i=2,Num_Images()
        integral = integral + integral[i]
      End Do
      ! Broadcast the result back to everyone.
      Do i=2,Num_Images()
        integral[i] = integral
      End Do
    End If
    Sync All
  End Subroutine
End Module