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


MODULE Erf_Series
   USE Precision

   IMPLICIT NONE
   PUBLIC
   INTEGER, PARAMETER :: max_iterations=100
   LOGICAL :: converged
   REAL(KIND=wp) :: error ! Error allowed in Erf(x) 

CONTAINS

   FUNCTION ErfReal(x) RESULT(erf)
       INTEGER :: status      ! I/O status
       INTEGER :: k ! Iteration counter
       REAL(KIND=wp) :: x, erf        ! x and erf(x)
       REAL(KIND=wp) :: sum, factor, expression     
       REAL(KIND=wp) :: epsilon    ! Error allowed in erf(x) 
       REAL(KIND=wp) :: pi     ! pi
   
      
          pi=4*ATAN(1.0_wp)    ! 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
                      converged=.TRUE.
                      erf=factor*sum
                      EXIT Summation
                    ELSE IF(k>=max_iterations) THEN        ! Not converged
                      converged=.FALSE.
                      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, erf is approximately 1
              erf=1.0_wp  
          END IF Magnitude      
    END FUNCTION ErfReal 
    
END MODULE Erf_Series    

PROGRAM Erf_Module
    USE Precision
    USE Erf_Series, ONLY: error, converged, ErfReal

    REAL(KIND=wp) :: x, erf       
    INTEGER :: status               ! i/o status - kicks you out if input is silly
 
    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
        
        erf=ErfReal(x)
        IF(converged) THEN
            WRITE(UNIT=*, FMT=*) "Erf(",x,") = ",erf
        ELSE
           WRITE(UNIT=*, FMT=*) " Convergence failed. Last estimate=", erf 
        END IF
   END DO MainLoop
END PROGRAM Erf_Module