function erfc03(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      real(wp)              :: t


      t = 1.0_wp - 7.5_wp/(abs(x)+3.75_wp)
      erfc = (((((((((((((((((((((+3.328130055126039e-10_wp *t    &
                                  -5.718639670776992e-10_wp)*t    &
                                  -4.066088879757269e-09_wp)*t    &
                                  +7.532536116142436e-09_wp)*t    &
                                  +3.026547320064576e-08_wp)*t    &
                                  -7.043998994397452e-08_wp)*t    &
                                  -1.822565715362025e-07_wp)*t    &
                                  +6.575825478226343e-07_wp)*t    &
                                  +7.478317101785790e-07_wp)*t    &
                                  -6.182369348098529e-06_wp)*t    &
                                  +3.584014089915968e-06_wp)*t    &
                                  +4.789838226695987e-05_wp)*t    &
                                  -1.524627476123466e-04_wp)*t    &
                                  -2.553523453642242e-05_wp)*t    &
                                  +1.802962431316418e-03_wp)*t    &
                                  -8.220621168415435e-03_wp)*t    &
                                  +2.414322397093253e-02_wp)*t    &
                                  -5.480232669380236e-02_wp)*t    &
                                  +1.026043120322792e-01_wp)*t    &
                                  -1.635718955239687e-01_wp)*t    &
                                  +2.260080669166197e-01_wp)*t    &
                                  -2.734219314954260e-01_wp)*t    &
                                  +1.455897212750385e-01_wp
      erfc = exp(-x**2)*erfc     
      if (x < 0.0_wp) erfc = 2.0_wp - erfc

      end function erfc03