diff --git a/camb/_compilers.py b/camb/_compilers.py
index 7bbbc8c7..6c8b372d 100644
--- a/camb/_compilers.py
+++ b/camb/_compilers.py
@@ -26,10 +26,10 @@ def get_ifort_version():
     return call_command("ifort -v")
 
 
-def get_gfortran_version():
-    ver = call_command("gfortran -dumpversion")
+def get_gfortran_version(command='gfortran'):
+    ver = call_command(command + " -dumpversion")
     if ver and '.' not in ver:
-        ver = call_command("gfortran -dumpfullversion")
+        ver = call_command(command + " -dumpfullversion")
     return ver
 
 
@@ -46,24 +46,34 @@ def check_gfortran(version=gfortran_min, msg=False, retry=False):
     else:
         ok = False
     if not ok and is_windows and not retry:
-        mingw = os.path.join(os.environ["ProgramFiles"], 'mingw-w64')
+        newpath = None
         best_version = gfortran_min
-        if os.path.isdir(mingw):
-            # look for mingw installation
-            dirs = [name for name in os.listdir(mingw) if
-                    gfortran_bits in name and os.path.isdir(os.path.join(mingw, name))]
-            path = None
-            for i, x in enumerate(dirs):
-                ver = x.split('-')[1]
-                if parse_version(best_version) <= parse_version(ver):
-                    best_version = ver
-                    path = x
-            if path:
-                newpath = os.path.join(mingw, path, 'mingw64', 'bin')
-                if os.path.exists(os.path.join(newpath, 'gfortran.exe')):
-                    if not compiler_environ["PATH"].startswith(newpath):
-                        compiler_environ["PATH"] = newpath + ';' + compiler_environ["PATH"]
-                    return check_gfortran(version, msg, retry=True)
+        for root in (os.path.expanduser('~'), os.environ["ProgramFiles"]):
+            mingw = os.path.join(root, 'mingw-w64')
+            if not os.path.isdir(mingw):
+                mingw = os.path.join(root, 'mingw64')
+
+            if os.path.isdir(mingw):
+                # look for mingw installation
+                dirs = [name for name in os.listdir(mingw) if
+                        gfortran_bits in name and os.path.isdir(os.path.join(mingw, name))]
+                for i, x in enumerate(dirs):
+                    if '.' in x:
+                        ver = x.split('-')[1]
+                        if parse_version(best_version) <= parse_version(ver):
+                            best_version = ver
+                            newpath = os.path.join(mingw, x, 'mingw64', 'bin')
+                bin = os.path.join(mingw, 'bin')
+                if os.path.exists(bin):
+                    ver = get_gfortran_version('"' + os.path.join(bin, 'gfortran') + '"')
+                    if ver and parse_version(best_version) <= parse_version(ver):
+                        best_version = ver
+                        newpath = bin
+        if newpath:
+            if os.path.exists(os.path.join(newpath, 'gfortran.exe')):
+                if not compiler_environ["PATH"].startswith(newpath):
+                    compiler_environ["PATH"] = newpath + ';' + compiler_environ["PATH"]
+                return check_gfortran(version, msg, retry=True)
     if ok and is_windows:
         version_str = str(subprocess.check_output("gfortran -dumpmachine", shell=True, env=compiler_environ))
         ok = gfortran_bits in version_str
diff --git a/camb/cambdll.dll b/camb/cambdll.dll
index 07a16182..29cabf8b 100644
Binary files a/camb/cambdll.dll and b/camb/cambdll.dll differ
diff --git a/fortran/DarkAge21cm.f90 b/fortran/DarkAge21cm.f90
index 1aee19e8..118be7ff 100644
--- a/fortran/DarkAge21cm.f90
+++ b/fortran/DarkAge21cm.f90
@@ -1,5 +1,6 @@
     module DarkAge21cm
     use constants
+    use splines
     !Functions for collision rates.
     implicit none
     contains
@@ -30,7 +31,7 @@ function kappa_HH_21cm(T, deriv)
         allocate(logRates(n_table),logTemps(n_table),ddlogRates(n_table))
         logRates = log(real(rates,dl)*0.01**3)
         logTemps = log(real(Temps,dl))
-        call spline(logTemps,logRates,n_table,1d30,1d30,ddlogRates)
+        call spline_def(logTemps,logRates,n_table,ddlogRates)
     end if
 
     if (T<=Temps(1)) then
@@ -135,7 +136,7 @@ function kappa_pH_21cm(T, deriv) ! from astro-ph/0702487
         allocate(logRates(n_table),logTemps(n_table),ddlogRates(n_table))
         logRates = log(real(rates,dl)*factor)
         logTemps = log(real(Temps,dl))
-        call spline(logTemps,logRates,n_table,1d30,1d30,ddlogRates)
+        call spline_def(logTemps,logRates,n_table,ddlogRates)
     end if
 
     if (T<=Temps(1)) then
diff --git a/fortran/DarkEnergyQuintessence.f90 b/fortran/DarkEnergyQuintessence.f90
index f7fcf222..631f3ee2 100644
--- a/fortran/DarkEnergyQuintessence.f90
+++ b/fortran/DarkEnergyQuintessence.f90
@@ -22,6 +22,7 @@ module Quintessence
     use results
     use constants
     use classes
+    use Interpolation
     implicit none
     private
 
diff --git a/fortran/SecondOrderPK.f90 b/fortran/SecondOrderPK.f90
index 17e3fa4a..0151e25f 100644
--- a/fortran/SecondOrderPK.f90
+++ b/fortran/SecondOrderPK.f90
@@ -176,6 +176,7 @@ end subroutine TSecondOrderPK_GetNonLinRatios_all
     subroutine GetRatios(this,CAMB_Pk, DoVel)
     !Fill the CAMB_Pk%nonlin_scaling array with sqrt(non-linear power/linear power)
     !for each redshift and wavenumber
+    use splines
     class(TSecondOrderPK) :: this
     type(MatterPowerData), target :: CAMB_Pk
     logical, intent(in):: DoVel
@@ -212,7 +213,7 @@ subroutine GetRatios(this,CAMB_Pk, DoVel)
 
             allocate(dPdLogK(CAMB_PK%num_k))
             !Get first derivative needed for series expansion at low r
-            call spline_deriv(CAMB_Pk%log_kh, CAMB_Pk%matpower(1,it), CAMB_Pk%ddmat(1,it), dPdLogK,CAMB_PK%num_k)
+            call spline_deriv(CAMB_Pk%log_kh, CAMB_Pk%matpower(:,it), CAMB_Pk%ddmat(:,it), dPdLogK,CAMB_PK%num_k)
 
             do i=1, CAMB_PK%num_k
 
diff --git a/fortran/SeparableBispectrum.f90 b/fortran/SeparableBispectrum.f90
index 090ea92c..41f1221b 100644
--- a/fortran/SeparableBispectrum.f90
+++ b/fortran/SeparableBispectrum.f90
@@ -20,6 +20,7 @@ module Bispectrum
     use SpherBessels
     use constants
     use MpiUtils
+    use splines
     implicit none
 
     integer, parameter :: max_bispectrum_deltas = 5, max_bispectrum_fields=3
@@ -78,7 +79,7 @@ subroutine InitBesselDerivs(CTrans)
     do i=1, CTrans%ls%nl
         !Spline agrees well
         !  call spline_deriv(BessRanges%points,ajl(1,i),ajlpr(1,i),dJl(1,i),BessRanges%npoints)
-        !  call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))
+        !  call spline_def(BessRanges%points,dJl(1,i),BessRanges%npoints,dddJl(1,i))
 
         l1 = CTrans%ls%l(i)
         do j=1, BessRanges%npoints
@@ -86,7 +87,7 @@ subroutine InitBesselDerivs(CTrans)
             call BJL(l1+1,BessRanges%points(j),Jp)
             dJl(j,i)=  ( l1*Jm - (l1+1)*Jp)/(2*l1+1)
         end do
-        call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))
+        call spline_def(BessRanges%points,dJl(:,i),BessRanges%npoints,dddJl(:,i))
 
     end do
 
diff --git a/fortran/bessels.f90 b/fortran/bessels.f90
index b88c8f01..e4704197 100644
--- a/fortran/bessels.f90
+++ b/fortran/bessels.f90
@@ -11,6 +11,7 @@ module SpherBessels
     use results
     use RangeUtils
     use MpiUtils
+    use splines
     implicit none
     private
 
@@ -115,8 +116,8 @@ subroutine GenerateBessels(lSamp, CP)
         end do
 
         !     get the interpolation matrix for bessel functions
-        call spline(BessRanges%points,ajl(1,j),num_xx,spl_large,spl_large,ajlpr(1,j))
-        call spline(BessRanges%points,ajlpr(1,j),num_xx,spl_large,spl_large,ddajlpr(1,j))
+        call spline_def(BessRanges%points,ajl(:,j),num_xx,ajlpr(:,j))
+        call spline_def(BessRanges%points,ajlpr(:,j),num_xx,ddajlpr(:,j))
     end do
     !$OMP END PARALLEL DO
 
diff --git a/fortran/camb_python.f90 b/fortran/camb_python.f90
index 272d6bbf..6d3db5ab 100644
--- a/fortran/camb_python.f90
+++ b/fortran/camb_python.f90
@@ -8,6 +8,7 @@ module handles
     use DarkEnergyPPF
     use ObjectLists
     use classes
+    use Interpolation
     implicit none
 
     Type c_MatterTransferData
@@ -284,10 +285,9 @@ end subroutine CAMBdata_GetSigma8
 
     subroutine CAMBdata_GetSigmaRArray(State, sigma, R, nR, z_ix, nz, var1, var2)
     Type(CAMBdata) :: State
+    integer, intent(in) :: nR, nz  
     real(dl) :: sigma(nR,nz)
-    integer :: nR
     real(dl) :: R(nR)
-    integer nz
     integer var1, var2, z_ix(nz)
     integer i, ix
 
@@ -303,8 +303,8 @@ end subroutine CAMBdata_GetSigmaRArray
 
     subroutine CAMBdata_GetMatterTransferks(State, nk, ks)
     Type(CAMBdata) :: State
-    real(dl) ks(nk)
     integer nk
+    real(dl) ks(nk)
 
     if (nk/=0) then
         ks =  State%MT%TransferData(Transfer_kh,:,1) * (State%CP%H0/100)
@@ -677,6 +677,7 @@ end function CAMB_TimeEvolution
 
 
     subroutine GetBackgroundThermalEvolution(this, ntimes, times, outputs)
+    use splines
     Type(CAMBdata) :: this
     integer, intent(in) :: ntimes
     real(dl), intent(in) :: times(ntimes)
diff --git a/fortran/cmbmain.f90 b/fortran/cmbmain.f90
index 0c91f233..fe3aafc1 100644
--- a/fortran/cmbmain.f90
+++ b/fortran/cmbmain.f90
@@ -754,7 +754,8 @@ subroutine GetSourceMem
         deallocate(ThisSources%LinearSrc)
     allocate(ThisSources%LinearSrc(ThisSources%Evolve_q%npoints,&
         ThisSources%SourceNum,State%TimeSteps%npoints), source=0._dl, stat=err)
-    if (err/=0) call GlobalError('Sources requires too much memory to allocate', error_unsupported_params)                                                                               
+    if (err/=0) call GlobalError('Sources requires too much memory to allocate', &
+        error_unsupported_params)                                                                               
 
     end subroutine GetSourceMem
 
@@ -1181,8 +1182,7 @@ subroutine MakeNonlinearSources
             !Do not use an associate for scaling. It does not work.
             scaling = State%CAMB_Pk%nonlin_ratio(ik,1:State%num_transfer_redshifts)
             if (all(abs(scaling-1) < 5e-4)) cycle
-            call spline(State%Transfer_Times(1), scaling(1), State%num_transfer_redshifts,&
-                spl_large, spl_large, ddScaling(1))
+            call spline_def(State%Transfer_Times, scaling, State%num_transfer_redshifts,ddScaling)
 
             tf_lo=1
             tf_hi=tf_lo+1
@@ -1222,8 +1222,8 @@ subroutine InitSourceInterpolation
     !$OMP PARALLEL DO DEFAULT(SHARED), SCHEDULE(STATIC), PRIVATE(i,j)
     do  i=1,State%TimeSteps%npoints
         do j=1, ThisSources%SourceNum
-            call spline(ThisSources%Evolve_q%points,ScaledSrc(1,j,i), &
-                ThisSources%Evolve_q%npoints,spl_large, spl_large, ddScaledSrc(1,j,i))
+            call spline_def(ThisSources%Evolve_q%points,ScaledSrc(:,j,i), &
+                ThisSources%Evolve_q%npoints, ddScaledSrc(:,j,i))
         end do
     end do
     !$OMP END PARALLEL DO
@@ -1374,8 +1374,8 @@ subroutine InterpolateSources(IV)
 
     if (.not.State%flat) then
         do i=1, ThisSources%SourceNum
-            call spline(State%TimeSteps%points,IV%Source_q(1,i),State%TimeSteps%npoints,&
-                spl_large,spl_large,IV%ddSource_q(1,i))
+            call spline_def(State%TimeSteps%points,IV%Source_q(:,i),State%TimeSteps%npoints,&
+                IV%ddSource_q(:,i))
         end do
     end if
 
diff --git a/fortran/config.f90 b/fortran/config.f90
index e1757fe8..bbf2d1f2 100644
--- a/fortran/config.f90
+++ b/fortran/config.f90
@@ -41,9 +41,6 @@ module config
 
     real(dl), parameter :: tol=1.0d-4 !Base tolerance for perturbation integrations
 
-    !     used as parameter for spline - tells it to use 'natural' end values
-    real(dl), parameter :: spl_large=1.e40_dl
-
     character(LEN=1024) :: highL_unlensed_cl_template = 'HighLExtrapTemplate_lenspotentialCls.dat'
     !fiducial high-accuracy high-L C_L used for making small cosmology-independent numerical corrections
     !to lensing and C_L interpolation. Ideally close to models of interest, but dependence is weak.
diff --git a/fortran/lensing.f90 b/fortran/lensing.f90
index e70be1ca..7f721cd5 100644
--- a/fortran/lensing.f90
+++ b/fortran/lensing.f90
@@ -39,6 +39,7 @@ module lensing
     use Precision
     use results
     use constants, only : const_pi, const_twopi, const_fourpi
+    use splines
     implicit none
     integer, parameter :: lensing_method_curv_corr=1,lensing_method_flat_corr=2, &
         lensing_method_harmonic=3
@@ -462,7 +463,7 @@ subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
             if (.false.) then
                 if (abs(sum(corrcontribs(1:jmax,1)))>1e-11) print *,i,sum(corrcontribs(1:jmax,1))
                 do j=1,4
-                    call spline(xl,corrcontribs(1,j),jmax,1d30,1d30,ddcontribs(1,j))
+                    call spline_def(xl,corrcontribs(:,j),jmax,ddcontribs(:,j))
                 end do
                 corr=0
                 llo=1
@@ -1157,7 +1158,7 @@ subroutine GetBessels(MaxArg)
             Bess(i,:)=Bessel_jn(0, maxbessel,x(i))
         end do
         do ix=0,maxbessel
-            call spline(x,Bess(1,ix),max_bes_ix,spl_large,spl_large,ddBess(1,ix))
+            call spline_def(x,Bess(:,ix),max_bes_ix,ddBess(:,ix))
         end do
 
         deallocate(x)
diff --git a/fortran/massive_neutrinos.f90 b/fortran/massive_neutrinos.f90
index 454d30a3..132ef913 100644
--- a/fortran/massive_neutrinos.f90
+++ b/fortran/massive_neutrinos.f90
@@ -108,6 +108,7 @@ subroutine TNuPerturbations_init(this,Accuracy)
     end subroutine TNuPerturbations_init
 
     subroutine ThermalNuBackground_init(this)
+    use splines
     class(TThermalNuBackground) :: this
     !  Initialize interpolation tables for massive neutrinos.
     !  Use cubic splines interpolation of log rhonu and pnu vs. log a*m.
@@ -141,6 +142,7 @@ end subroutine ThermalNuBackground_init
     subroutine nuRhoPres(am,rhonu,pnu)
     !  Compute the density and pressure of one eigenstate of massive neutrinos,
     !  in units of the mean density of one flavor of massless neutrinos.
+    use splines
     real(dl),  parameter :: qmax=30._dl
     integer, parameter :: nq=100
     real(dl) dum1(nq+1),dum2(nq+1)
diff --git a/fortran/recfast.f90 b/fortran/recfast.f90
index d72836b9..4ff52308 100644
--- a/fortran/recfast.f90
+++ b/fortran/recfast.f90
@@ -501,7 +501,7 @@ subroutine TRecfast_init(this,State, WantTSpin)
     implicit none
     class(TRecfast), target :: this
     class(TCAMBdata), target :: State
-    real(dl) :: Trad,Tmat,Tspin,d0hi,d0lo
+    real(dl) :: Trad,Tmat,Tspin
     integer :: I
     Type(RecombinationData), pointer :: Calc
     logical, intent(in), optional :: WantTSpin
@@ -701,12 +701,10 @@ subroutine TRecfast_init(this,State, WantTSpin)
             !          write (*,'(5E15.5)') zend, Trad, Tmat, Tspin, x
         end do
 
-        d0hi=1.0d40
-        d0lo=1.0d40
-        call spline(Calc%zrec,Calc%xrec,nz,d0lo,d0hi,Calc%dxrec)
-        call spline(Calc%zrec,Calc%tmrec,nz,d0lo,d0hi,Calc%dtmrec)
+        call spline_def(Calc%zrec,Calc%xrec,nz,Calc%dxrec)
+        call spline_def(Calc%zrec,Calc%tmrec,nz,Calc%dtmrec)
         if (Calc%doTspin) then
-            call spline(Calc%zrec,Calc%tsrec,nz,d0lo,d0hi,Calc%dtsrec)
+            call spline_def(Calc%zrec,Calc%tsrec,nz,Calc%dtsrec)
         end if
     class default
         call MpiStop('Wrong state type')
diff --git a/fortran/results.f90 b/fortran/results.f90
index a120100a..9c59dc91 100644
--- a/fortran/results.f90
+++ b/fortran/results.f90
@@ -31,6 +31,7 @@ module results
     use MathUtils
     use config
     use model
+    use splines
     implicit none
     public
 
@@ -96,7 +97,7 @@ module results
 
     !Sources
     Type CalWins
-        real(dl), allocatable :: awin_lens(:),  dawin_lens(:)
+        real(dl), allocatable :: awin_lens(:), dawin_lens(:)
     end Type CalWins
 
     Type LimberRec
@@ -606,7 +607,7 @@ function CAMBdata_DeltaTime(this, a1,a2, in_tol)
 
     end function CAMBdata_DeltaTime
 
-    subroutine CAMBdata_DeltaTimeArr(this, arr,  a1, a2, n, tol)
+    subroutine CAMBdata_DeltaTimeArr(this, arr, a1, a2, n, tol)
     class(CAMBdata) :: this
     integer, intent(in) :: n
     real(dl), intent(out) :: arr(n)
@@ -630,7 +631,7 @@ function CAMBdata_TimeOfz(this, z, tol)
     CAMBdata_TimeOfz= this%DeltaTime(0._dl,1._dl/(z+1._dl), tol)
     end function CAMBdata_TimeOfz
 
-    subroutine CAMBdata_TimeOfzArr(this, arr,  z, n, tol)
+    subroutine CAMBdata_TimeOfzArr(this, arr, z, n, tol)
     !z array must be monotonically *decreasing* so times increasing
     class(CAMBdata) :: this
     integer, intent(in) :: n
@@ -671,7 +672,7 @@ function CAMBdata_DeltaPhysicalTimeGyr(this, a1,a2, in_tol)
     CAMBdata_DeltaPhysicalTimeGyr = Integrate_Romberg(this, dtda,a1,a2,atol)*Mpc/c/Gyr
     end function CAMBdata_DeltaPhysicalTimeGyr
 
-    subroutine CAMBdata_DeltaPhysicalTimeGyrArr(this, arr,  a1, a2, n, tol)
+    subroutine CAMBdata_DeltaPhysicalTimeGyrArr(this, arr, a1, a2, n, tol)
     class(CAMBdata) :: this
     integer, intent(in) :: n
     real(dl), intent(out) :: arr(n)
@@ -1499,14 +1500,13 @@ end subroutine lSamples_init
 
     subroutine InterpolateClArr(lSet,iCl, all_Cl, max_index)
     class(lSamples), intent(in) :: lSet
-    real(dl), intent(in) :: iCl(*)
+    real(dl), intent(in) :: iCl(1:*)
     real(dl), intent(out):: all_Cl(lSet%lmin:*)
     integer, intent(in), optional :: max_index
     integer il,llo,lhi, xi
     real(dl) ddCl(lSet%nl)
     real(dl) xl(lSet%nl)
     real(dl) a0,b0,ho
-    real(dl), parameter :: cllo=1.e30_dl,clhi=1.e30_dl
     integer max_ind
 
     max_ind = PresentDefault(lSet%nl, max_index)
@@ -1514,7 +1514,7 @@ subroutine InterpolateClArr(lSet,iCl, all_Cl, max_index)
     if (max_ind > lSet%nl) call MpiStop('Wrong max_ind in InterpolateClArr')
 
     xl = real(lSet%l(1:lSet%nl),dl)
-    call spline(xl,iCL(1),max_ind,cllo,clhi,ddCl(1))
+    call spline_def(xl,iCL,max_ind,ddCl)
 
     llo=1
     do il=lSet%lmin,lSet%l(max_ind)
@@ -2328,7 +2328,7 @@ subroutine SetTimeSteps(this,State,TimeSteps)
 
                 if (Win%kind /= window_lensing  .and. &
                     Win%tau_peakend-Win%tau_peakstart < nint(60*TimeSampleBoost) * delta) then
-                    call TimeSteps%Add(Win%tau_peakstart,Win%tau_peakend,  nint(60*TimeSampleBoost))
+                    call TimeSteps%Add(Win%tau_peakstart,Win%tau_peakend, nint(60*TimeSampleBoost))
                     !This should be over peak
                 end if
             end associate
@@ -2388,7 +2388,7 @@ subroutine SetTimeStepWindows(this,State,TimeSteps)
     class(CAMBdata) :: State
     Type(TRanges) :: TimeSteps
     integer i, j, jstart, ix
-    real(dl) tau,  a, a2
+    real(dl) tau, a, a2
     real(dl) Tspin, Trad, rho_fac, tau_eps
     real(dl) window, winamp
     real(dl) z,rhos, adot, exp_fac
@@ -2507,38 +2507,39 @@ subroutine SetTimeStepWindows(this,State,TimeSteps)
         associate (RedWin => State%Redshift_W(i))
 
             ! int (a*rho_s/H)' a W_f(a) d\eta, or for counts int g/chi deta
-            call spline(TimeSteps%points(jstart),int_tmp(jstart,i),ninterp,spl_large,spl_large,tmp)
-            call spline_integrate(TimeSteps%points(jstart),int_tmp(jstart,i),tmp, tmp2(jstart),ninterp)
+            call spline_def(TimeSteps%points(jstart:),int_tmp(jstart:,i),ninterp,tmp)
+            call spline_integrate(TimeSteps%points(jstart:),int_tmp(jstart:,i),tmp, tmp2(jstart:),ninterp)
             RedWin%WinV(jstart:TimeSteps%npoints) =  &
                 RedWin%WinV(jstart:TimeSteps%npoints) + tmp2(jstart:TimeSteps%npoints)
 
-            call spline(TimeSteps%points(jstart),RedWin%WinV(jstart),ninterp,spl_large,spl_large,RedWin%ddWinV(jstart))
-            call spline_deriv(TimeSteps%points(jstart),RedWin%WinV(jstart),RedWin%ddWinV(jstart), RedWin%dWinV(jstart), ninterp)
+            call spline_def(TimeSteps%points(jstart:),RedWin%WinV(jstart:),ninterp,RedWin%ddWinV(jstart:))
+            call spline_deriv(TimeSteps%points(jstart:),RedWin%WinV(jstart:),RedWin%ddWinV(jstart:), RedWin%dWinV(jstart:), ninterp)
 
-            call spline(TimeSteps%points(jstart),RedWin%Wing(jstart),ninterp,spl_large,spl_large,RedWin%ddWing(jstart))
-            call spline_deriv(TimeSteps%points(jstart),RedWin%Wing(jstart),RedWin%ddWing(jstart), RedWin%dWing(jstart), ninterp)
+            call spline_def(TimeSteps%points(jstart:),RedWin%Wing(jstart:),ninterp,RedWin%ddWing(jstart:))
+            call spline_deriv(TimeSteps%points(jstart:),RedWin%Wing(jstart:),RedWin%ddWing(jstart:), RedWin%dWing(jstart:), ninterp)
 
-            call spline(TimeSteps%points(jstart),RedWin%Wing2(jstart),ninterp,spl_large,spl_large,RedWin%ddWing2(jstart))
-            call spline_deriv(TimeSteps%points(jstart),RedWin%Wing2(jstart),RedWin%ddWing2(jstart), &
-                RedWin%dWing2(jstart), ninterp)
+            call spline_def(TimeSteps%points(jstart:),RedWin%Wing2(jstart:),ninterp,RedWin%ddWing2(jstart:))
+            call spline_deriv(TimeSteps%points(jstart:),RedWin%Wing2(jstart:),RedWin%ddWing2(jstart:), &
+                RedWin%dWing2(jstart:), ninterp)
 
-            call spline_integrate(TimeSteps%points(jstart),RedWin%Wing(jstart),RedWin%ddWing(jstart), RedWin%WinF(jstart),ninterp)
+            call spline_integrate(TimeSteps%points(jstart:),RedWin%Wing(jstart:), &
+                RedWin%ddWing(jstart:), RedWin%WinF(jstart:),ninterp)
             RedWin%Fq = RedWin%WinF(TimeSteps%npoints)
 
             if (RedWin%kind == window_21cm) then
-                call spline_integrate(TimeSteps%points(jstart),RedWin%Wing2(jstart),&
-                    RedWin%ddWing2(jstart), tmp(jstart),ninterp)
+                call spline_integrate(TimeSteps%points(jstart:),RedWin%Wing2(jstart:),&
+                    RedWin%ddWing2(jstart:), tmp(jstart:),ninterp)
                 RedWin%optical_depth_21 = tmp(TimeSteps%npoints) / (State%CP%TCMB*1000)
                 !WinF not used.. replaced below
 
-                call spline(TimeSteps%points(jstart),RedWin%Wingtau(jstart),ninterp,spl_large,spl_large,RedWin%ddWingtau(jstart))
-                call spline_deriv(TimeSteps%points(jstart),RedWin%Wingtau(jstart),RedWin%ddWingtau(jstart), &
-                    RedWin%dWingtau(jstart), ninterp)
+                call spline_def(TimeSteps%points(jstart:),RedWin%Wingtau(jstart:),ninterp,RedWin%ddWingtau(jstart:))
+                call spline_deriv(TimeSteps%points(jstart:),RedWin%Wingtau(jstart:),RedWin%ddWingtau(jstart:), &
+                    RedWin%dWingtau(jstart:), ninterp)
             elseif (RedWin%kind == window_counts) then
 
                 if (State%CP%SourceTerms%counts_evolve) then
-                    call spline(TimeSteps%points(jstart),back_count_tmp(jstart,i),ninterp,spl_large,spl_large,tmp)
-                    call spline_deriv(TimeSteps%points(jstart),back_count_tmp(jstart,i),tmp,tmp2(jstart),ninterp)
+                    call spline_def(TimeSteps%points(jstart:),back_count_tmp(jstart:,i),ninterp,tmp)
+                    call spline_deriv(TimeSteps%points(jstart:),back_count_tmp(jstart:,i),tmp,tmp2(jstart:),ninterp)
                     do ix = jstart, TimeSteps%npoints
                         if (RedWin%Wing(ix)==0._dl) then
                             RedWin%Wingtau(ix) = 0
@@ -2552,8 +2553,8 @@ subroutine SetTimeStepWindows(this,State,TimeSteps)
                     end do
 
                     !comoving_density_ev is d log(a^3 n_s)/d eta * window
-                    call spline(TimeSteps%points(jstart),RedWin%comoving_density_ev(jstart),ninterp,spl_large,spl_large,tmp)
-                    call spline_deriv(TimeSteps%points(jstart),RedWin%comoving_density_ev(jstart),tmp,tmp2(jstart),ninterp)
+                    call spline_def(TimeSteps%points(jstart:),RedWin%comoving_density_ev(jstart:),ninterp,tmp)
+                    call spline_deriv(TimeSteps%points(jstart:),RedWin%comoving_density_ev(jstart:),tmp,tmp2(jstart:),ninterp)
                     do ix = jstart, TimeSteps%npoints
                         if (RedWin%Wing(ix)==0._dl) then
                             RedWin%comoving_density_ev(ix) = 0
@@ -2565,8 +2566,8 @@ subroutine SetTimeStepWindows(this,State,TimeSteps)
                     end do
                 else
                     RedWin%comoving_density_ev=0
-                    call spline(TimeSteps%points(jstart),hubble_tmp(jstart),ninterp,spl_large,spl_large,tmp)
-                    call spline_deriv(TimeSteps%points(jstart),hubble_tmp(jstart),tmp, tmp2(jstart), ninterp)
+                    call spline_def(TimeSteps%points(jstart:),hubble_tmp(jstart:),ninterp,tmp)
+                    call spline_deriv(TimeSteps%points(jstart:),hubble_tmp(jstart:),tmp, tmp2(jstart:), ninterp)
 
                     !assume d( a^3 n_s) of background population is zero, so remaining terms are
                     !wingtau =  g*(2/H\chi + Hdot/H^2)  when s=0; int_tmp = window/chi
@@ -2578,14 +2579,14 @@ subroutine SetTimeStepWindows(this,State,TimeSteps)
                         *RedWin%Wing(jstart:TimeSteps%npoints)
                 endif
 
-                call spline(TimeSteps%points(jstart),RedWin%Wingtau(jstart),ninterp, &
-                    spl_large,spl_large,RedWin%ddWingtau(jstart))
-                call spline_deriv(TimeSteps%points(jstart),RedWin%Wingtau(jstart),RedWin%ddWingtau(jstart), &
-                    RedWin%dWingtau(jstart), ninterp)
+                call spline_def(TimeSteps%points(jstart:),RedWin%Wingtau(jstart:),ninterp, &
+                    RedWin%ddWingtau(jstart:))
+                call spline_deriv(TimeSteps%points(jstart:),RedWin%Wingtau(jstart:),RedWin%ddWingtau(jstart:), &
+                    RedWin%dWingtau(jstart:), ninterp)
 
                 !WinF is int[ g*(...)]
-                call spline_integrate(TimeSteps%points(jstart),RedWin%Wingtau(jstart),&
-                    RedWin%ddWingtau(jstart), RedWin%WinF(jstart),ninterp)
+                call spline_integrate(TimeSteps%points(jstart:),RedWin%Wingtau(jstart:),&
+                    RedWin%ddWingtau(jstart:), RedWin%WinF(jstart:),ninterp)
             end if
         end associate
     end do
@@ -2721,7 +2722,7 @@ subroutine Init_ClTransfer(CTrans)
     call CTrans%q%getArray(.true.)
 
     allocate(CTrans%Delta_p_l_k(CTrans%NumSources,&
-        min(CTrans%max_index_nonlimber,CTrans%ls%nl), CTrans%q%npoints),  STAT = st)
+        min(CTrans%max_index_nonlimber,CTrans%ls%nl), CTrans%q%npoints), STAT = st)
     if (st /= 0) call MpiStop('Init_ClTransfer: Error allocating memory for transfer functions')
     CTrans%Delta_p_l_k = 0
 
@@ -2790,7 +2791,7 @@ subroutine TCLdata_InitCls(this, State)
             allocate(this%Cl_scalar(CP%Min_l:CP%Max_l, C_Temp:State%Scalar_C_last), source=0._dl)
             if (CP%want_cl_2D_array) then
                 if (allocated(this%Cl_scalar_array)) deallocate(this%Cl_scalar_array)
-                allocate(this%Cl_scalar_Array(CP%Min_l:CP%Max_l,  &
+                allocate(this%Cl_scalar_Array(CP%Min_l:CP%Max_l, &
                     3+State%num_redshiftwindows+CP%CustomSources%num_custom_sources, &
                     3+State%num_redshiftwindows+CP%CustomSources%num_custom_sources))
                 this%Cl_scalar_array = 0
@@ -3013,14 +3014,14 @@ subroutine TCLdata_NormalizeClsAtL(this,CP,lnorm)
 
     if (CP%WantScalars) then
         Norm=1/this%Cl_scalar(lnorm, C_Temp)
-        this%Cl_scalar(CP%Min_l:CP%Max_l,  C_Temp:C_Cross) = this%Cl_scalar(CP%Min_l:CP%Max_l, C_Temp:C_Cross) * Norm
+        this%Cl_scalar(CP%Min_l:CP%Max_l, C_Temp:C_Cross) = this%Cl_scalar(CP%Min_l:CP%Max_l, C_Temp:C_Cross) * Norm
     end if
 
     if (CP%WantTensors) then
         if (.not.CP%WantScalars) Norm = 1/this%Cl_tensor(lnorm, C_Temp)
         !Otherwise Norm already set correctly
         this%Cl_tensor(CP%Min_l:CP%Max_l_tensor, CT_Temp:CT_Cross) =  &
-            this%Cl_tensor(CP%Min_l:CP%Max_l_tensor,  CT_Temp:CT_Cross) * Norm
+            this%Cl_tensor(CP%Min_l:CP%Max_l_tensor, CT_Temp:CT_Cross) * Norm
     end if
 
     end subroutine TCLdata_NormalizeClsAtL
@@ -3262,11 +3263,10 @@ end subroutine MatterPowerData_Load
     subroutine MatterPowerdata_getsplines(PK_data)
     Type(MatterPowerData) :: PK_data
     integer i
-    real(dl), parameter :: cllo=1.e30_dl,clhi=1.e30_dl
 
     do i = 1,PK_Data%num_z
-        call spline(PK_data%log_kh,PK_data%matpower(1,i),PK_data%num_k,&
-            cllo,clhi,PK_data%ddmat(1,i))
+        call spline_def(PK_data%log_kh,PK_data%matpower(:,i),PK_data%num_k,&
+            PK_data%ddmat(:,i))
     end do
 
     end subroutine MatterPowerdata_getsplines
@@ -3275,15 +3275,14 @@ end subroutine MatterPowerdata_getsplines
     subroutine MatterPowerdata_getsplines21cm(PK_data)
     Type(MatterPowerData) :: PK_data
     integer i
-    real(dl), parameter :: cllo=1.e30_dl,clhi=1.e30_dl
 
     do i = 1,PK_Data%num_z
-        call spline(PK_data%log_k,PK_data%matpower(1,i),PK_data%num_k,&
-            cllo,clhi,PK_data%ddmat(1,i))
-        call spline(PK_data%log_k,PK_data%vvpower(1,i),PK_data%num_k,&
-            cllo,clhi,PK_data%ddvvpower(1,i))
-        call spline(PK_data%log_k,PK_data%vdpower(1,i),PK_data%num_k,&
-            cllo,clhi,PK_data%ddvdpower(1,i))
+        call spline_def(PK_data%log_k,PK_data%matpower(:,i),PK_data%num_k,&
+            PK_data%ddmat(:,i))
+        call spline_def(PK_data%log_k,PK_data%vvpower(:,i),PK_data%num_k,&
+            PK_data%ddvvpower(:,i))
+        call spline_def(PK_data%log_k,PK_data%vdpower(:,i),PK_data%num_k,&
+            PK_data%ddvdpower(:,i))
     end do
 
     end subroutine MatterPowerdata_getsplines21cm
@@ -3308,7 +3307,7 @@ subroutine MatterPowerdata_Free(PK_data)
 
     end subroutine MatterPowerdata_Free
 
-    function MatterPowerData_k(PK,  kh, itf, index_cache) result(outpower)
+    function MatterPowerData_k(PK, kh, itf, index_cache) result(outpower)
     !Get matter power spectrum at particular k/h by interpolation
     Type(MatterPowerData) :: PK
     integer, intent(in) :: itf
@@ -3359,7 +3358,7 @@ function MatterPowerData_k(PK,  kh, itf, index_cache) result(outpower)
     end function MatterPowerData_k
 
     !Sources
-    subroutine MatterPower21cm_k(PK,  k, itf, monopole, vv, vd)
+    subroutine MatterPower21cm_k(PK, k, itf, monopole, vv, vd)
     !Get monopole and velocity power at particular k by interpolation
     Type(MatterPowerData) :: PK
     integer, intent(in) :: itf
@@ -3429,7 +3428,7 @@ subroutine Transfer_GetMatterPowerS(State, MTrans, outpower, itf, minkh, dlnkh,
     real(dl):: minkhd, dlnkhd
 
     minkhd = minkh; dlnkhd = dlnkh
-    call Transfer_GetMatterPowerD(State, MTrans,  outpowerd, itf, minkhd, dlnkhd, npoints,var1, var2)
+    call Transfer_GetMatterPowerD(State, MTrans, outpowerd, itf, minkhd, dlnkhd, npoints,var1, var2)
     outpower(1:npoints) = outpowerd(1:npoints)
 
     end subroutine Transfer_GetMatterPowerS
@@ -3456,7 +3455,6 @@ subroutine Transfer_GetMatterPowerD(State, MTrans, outpower, itf_PK, minkh, dlnk
     real(dl), intent(in) :: minkh, dlnkh
     integer, intent(in), optional :: var1, var2
 
-    real(dl), parameter :: cllo=1.e30_dl,clhi=1.e30_dl
     integer ik, llo,il,lhi,lastix
     real(dl) matpower(MTrans%num_q_trans), kh, kvals(MTrans%num_q_trans), ddmat(MTrans%num_q_trans)
     real(dl) atransfer,xi, a0, b0, ho, logmink,k, h
@@ -3503,7 +3501,7 @@ subroutine Transfer_GetMatterPowerD(State, MTrans, outpower, itf_PK, minkh, dlnk
         endif
     endif
     if (log_interp) matpower = log(sign*matpower)
-    call spline(kvals,matpower,MTrans%num_q_trans,cllo,clhi,ddmat)
+    call spline_def(kvals,matpower,MTrans%num_q_trans, ddmat)
 
     llo=1
     lastix = npoints + 1
@@ -3712,7 +3710,7 @@ subroutine Transfer_Get_sigma8(State, MTrans, R, var1, var2)
 
     end subroutine Transfer_Get_sigma8
 
-    subroutine Transfer_Get_sigmas(State, MTrans,  R, var_delta, var_v)
+    subroutine Transfer_Get_sigmas(State, MTrans, R, var_delta, var_v)
     !Get sigma8 and sigma_{delta v} (for growth, like f sigma8 in LCDM)
     class(CAMBdata) :: State
     Type(MatterTransferData) :: MTrans
@@ -3727,7 +3725,7 @@ subroutine Transfer_Get_sigmas(State, MTrans,  R, var_delta, var_v)
     s1 = PresentDefault (transfer_power_var, var_delta)
     s2 = PresentDefault (Transfer_Newt_vel_cdm, var_v)
 
-    call Transfer_Get_SigmaR(State, MTrans,  radius, MTrans%sigma_8, s1,s1)
+    call Transfer_Get_SigmaR(State, MTrans, radius, MTrans%sigma_8, s1,s1)
     if (State%get_growth_sigma8) call Transfer_Get_SigmaR(State, MTrans, radius, &
         MTrans%sigma2_vdelta_8(:), s1, s2, root=.false.)
 
@@ -3864,7 +3862,7 @@ subroutine Transfer_SaveMatterPower(MTrans, State,FileNames, all21cm)
                 points = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/dlnkh+1
                 !             dlnkh = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/(points-0.999)
                 allocate(outpower(points,1))
-                call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)
+                call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf, minkh,dlnkh, points)
 
                 columns(1) = 'P'
                 unit = open_file_header(FileNames(itf), 'k/h', columns(:1), 15)
@@ -3956,7 +3954,7 @@ function Get21cmCl_l(Vars,kin)
     end if
     x= Vars%chi*k
 
-    call MatterPower21cm_k(Vars%PK,  k, Vars%itf, monopole, vv, vd)
+    call MatterPower21cm_k(Vars%PK, k, Vars%itf, monopole, vv, vd)
     call bjl_external(Vars%l, x, jl)
     call bjl_external(Vars%l-1, x, jlm1)
     ddjl = -( 2/x*jlm1-(Vars%l+2)*real(Vars%l+1,dl)/x**2*jl + jl)
@@ -3982,7 +3980,7 @@ function Get21cmCl_l_avg(Vars,kin)
     end if
     x= Vars%chi*k
 
-    call MatterPower21cm_k(Vars%PK,  k, Vars%itf, monopole, vv, vd)
+    call MatterPower21cm_k(Vars%PK, k, Vars%itf, monopole, vv, vd)
     lphalf=Vars%l+0.5_dl
 
     jl = 1/(2*x**2) /sqrt(1-(lphalf/x)**2)
diff --git a/fortran/subroutines.f90 b/fortran/subroutines.f90
index 4974ac6e..fdca4d5c 100644
--- a/fortran/subroutines.f90
+++ b/fortran/subroutines.f90
@@ -1,11 +1,25 @@
     !Low-level numerical routines for splines and dverk for differential equation integration.
 
-    subroutine splder(y,dy,n, g)
+    module splines
     use Precision
+    use Interpolation
+    implicit none
+    contains
+    
+    subroutine spline_def(x,y,n,d2)
+    !Low-level initialize spline arrays with default boundary conditions 
+    integer, intent(in) :: n
+    real(sp_acc), intent(in) :: x(n), y(n)
+    real(sp_acc), intent(out) :: d2(n)
+    
+    call spline(x,y,n,SPLINE_DANGLE,SPLINE_DANGLE,d2)
+    
+    end subroutine spline_def
+
+    subroutine splder(y,dy,n, g)
     !  Splder fits a cubic spline to y and returns the first derivatives at
     !  the grid points in dy.  Dy is equivalent to a 4th-order Pade
     !  difference formula for dy/di.
-    implicit none
     integer, intent(in) :: n
     real(dl), intent(in) :: y(n),g(n)
     real(dl), intent(out) :: dy(n)
@@ -31,9 +45,7 @@ subroutine splder(y,dy,n, g)
     end subroutine splder
     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     subroutine splini(g,n)
-    use Precision
     !  Splini must be called before splder to initialize array g in common.
-    implicit none
     integer, intent(in) :: n
     real(dl), intent(out):: g(n)
     integer :: i
@@ -44,61 +56,8 @@ subroutine splini(g,n)
     end do
     end subroutine splini
 
-
-    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-    ! calculates array of second derivatives used by cubic spline
-    ! interpolation. y2 is array of second derivatives, yp1 and ypn are first
-    ! derivatives at end points.
-
-    !Thanks Martin Reinecke
-    subroutine spline(x,y,n,d11,d1n,d2)
-    use Precision
-    integer, intent(in) :: n
-    real(dl), intent(in) :: x(n), y(n), d11, d1n
-    real(dl), intent(out) :: d2(n)
-    integer i
-    real(dl) xp,qn,sig,un,xxdiv,u(n-1),d1l,d1r
-
-    d1r= (y(2)-y(1))/(x(2)-x(1))
-    if (d11>.99e30_dl) then
-        d2(1)=0._dl
-        u(1)=0._dl
-    else
-        d2(1)=-0.5_dl
-        u(1)=(3._dl/(x(2)-x(1)))*(d1r-d11)
-    endif
-
-    do i=2,n-1
-        d1l=d1r
-        d1r=(y(i+1)-y(i))/(x(i+1)-x(i))
-        xxdiv=1._dl/(x(i+1)-x(i-1))
-        sig=(x(i)-x(i-1))*xxdiv
-        xp=1._dl/(sig*d2(i-1)+2._dl)
-
-        d2(i)=(sig-1._dl)*xp
-
-        u(i)=(6._dl*(d1r-d1l)*xxdiv-sig*u(i-1))*xp
-    end do
-    d1l=d1r
-
-    if (d1n>.99e30_dl) then
-        qn=0._dl
-        un=0._dl
-    else
-        qn=0.5_dl
-        un=(3._dl/(x(n)-x(n-1)))*(d1n-d1l)
-    endif
-
-    d2(n)=(un-qn*u(n-1))/(qn*d2(n-1)+1._dl)
-    do i=n-1,1,-1
-        d2(i)=d2(i)*d2(i+1)+u(i)
-    end do
-    end subroutine spline
-
     SUBROUTINE spline_deriv(x,y,y2,y1,n)
     !Get derivative y1 given array of x, y and y''
-    use Precision
-    implicit none
     INTEGER, intent(in) :: n
     real(dl), intent(in) :: x(n), y(n), y2(n)
     real(dl), intent(out) :: y1(n)
@@ -106,7 +65,6 @@ SUBROUTINE spline_deriv(x,y,y2,y1,n)
     real(dl) dx
 
     do i=1, n-1
-
         dx = (x(i+1) - x(i))
         y1(i) = (y(i+1) - y(i))/dx - dx*(2*y2(i) + y2(i+1))/6
     end do
@@ -117,7 +75,6 @@ END SUBROUTINE spline_deriv
 
     subroutine spline_integrate(x,y,y2,yint,n)
     !Cumulative integral of cubic spline
-    use Precision
     integer, intent(in) :: n
     real(dl), intent(in) :: x(n), y(n), y2(n)
     real(dl), intent(out) :: yint(n)
@@ -140,12 +97,10 @@ end subroutine spline_integrate
 
 
     subroutine splint(y,z,n)
-    use Precision
     !  Splint integrates a cubic spline, providing the output value
     !  z = integral from 1 to n of s(i)di, where s(i) is the spline fit
     !  to y(i).
     !
-    implicit none
     integer, intent(in) :: n
     real(dl), intent(in) :: y(n)
     real(dl), intent(out) :: z
@@ -163,6 +118,8 @@ subroutine splint(y,z,n)
     z= z + sum(y(2:n1))
     end subroutine splint
 
+    end module splines
+
 
     !This version is modified to pass an object parameter to the function on each call
     !Fortunately Fortran doesn't do type checking on functions, so we can pretend the
@@ -925,6 +882,6 @@ subroutine dverk (EV,n, fcn, x, y, xend, tol, ind, c, nw, w)
     !
     end subroutine dverk
 
-  
+
 
 
diff --git a/setup.py b/setup.py
index adc60489..abb10777 100644
--- a/setup.py
+++ b/setup.py
@@ -124,6 +124,7 @@ def make_library(cluster=False):
                   'you may need to log off and on again).' % _compile.gfortran_min)
             print('        You can get a Windows gfortran build from https://sourceforge.net/projects/mingw-w64/files/')
             print('        - go to Files, and download MinGW-W64 Online Installer.')
+            print('        Alternatively news versions at https://github.com/niXman/mingw-builds-binaries')
             if _compile.is_32_bit:
                 raise IOError('No 32bit Windows DLL provided, you need to build or use 64 bit python')
             else: