MODULE Precision
   PUBLIC
   INTEGER, PARAMETER :: sp=KIND(0.0E0), dp=KIND(0.0D0)
   INTEGER, PARAMETER :: wp=sp ! Use single precision in this run
END MODULE Precision


PROGRAM Erf_Series
   USE Precision

   IMPLICIT NONE
   INTEGER, PARAMETER :: max_iterations=100
   INTEGER :: status ! I/O status
   INTEGER :: k ! Iteration counter
   REAL(KIND=wp) :: x, erf ! x and erf(x)
   REAL(KIND=wp) :: sum, factor, expression ! A few auxiliary variables
   REAL(KIND=wp) :: error, epsilon ! Error allowed in erf(x) and in the sum
   REAL(KIND=wp) :: pi ! The constant pi
   
   MainLoop: DO
   
      WRITE(UNIT=*,FMT="(A)",ADVANCE="NO") "Enter the value of x and epsilon ('S' to stop): "
      READ(UNIT=*,FMT=*,IOSTAT=status) x,error      
      IF(status/=0) EXIT MainLoop

      pi=4*ATAN(1.0_wp) ! This is a good way of calculating pi to the needed precision
      factor=2*x/SQRT(pi) ! The factor in front of the sum
      epsilon=error/factor ! This is the error allowed *in* the sum     

      Magnitude: IF(ABS(x)<5.0_wp) THEN ! We can use the series provided
        sum=1.0_wp ! This is E(0)
        expression=1.0_wp ! This is E(0) as well 
        k=0 ! Initialize counter
        Summation: DO
          k=k+1 ! Iteration counter
          ! Calculate the next term to be added
          expression=-REAL(2*k-1,wp)/REAL(2*k+1,wp)/REAL(k,wp)*(x**2)*expression
          ! Now check for convergence:
          IF(ABS(expression)<epsilon) THEN ! Converged
            WRITE(UNIT=*,FMT=*) "Convergence was achieved with ",k," terms in the sum"
            erf=factor*sum
            WRITE(UNIT=*,FMT=*) "Erf (",x,") = ",erf
            EXIT Summation
          ELSE IF(k>=max_iterations) THEN ! Not converged
            WRITE(UNIT=*,FMT=*) "Convergence was not achieved!"   
            EXIT Summation
          ELSE ! Keep on adding terms
            sum=sum+expression ! Add the term to the sum
            CYCLE Summation ! This can be ommitted
          END IF
        END DO Summation
      ELSE ! If x is large in magnitude, we can not use the above series easily
        erf=1.0_wp
        WRITE(UNIT=*,FMT=*) "The value of x is too large! An approximate answer is:"
        WRITE(UNIT=*,FMT=*) "Erf (",x,") = ",erf     
      END IF Magnitude    
   END DO MainLoop
   
END PROGRAM Erf_Series