better time step adaption for sharp features in windows

camb/sources.py

@@ -6,6 +6,8 @@ class SourceWindow(F2003Class):
     """
     Abstract base class for a number count/lensing/21cm source window function.
     A list of instances of these classes can be assigned to the SourceWindows field of :class:`.model.CAMBparams`.
+
+    Note that source windows can currently only be used in flat models.
     """
     _fields_ = [("source_type", c_int, {"names": ["21cm", "counts", "lensing"], "start": 1}),
                 ("bias", c_double),

fortran/SeparableBispectrum.f90

@@ -366,7 +366,6 @@ subroutine GetBispectrum(State,CTrans)
             end if
         end if
         SampleL%nl= SampleL%nl + 1
-        ! print *,l1
         SampleL%l(SampleL%nl) = l1
         if (l1 == lmax) exit
     end do

fortran/VisualStudio/CAMBLib.vfproj

@@ -31,14 +31,12 @@
 		<File RelativePath="..\cmbmain.f90"/>
 		<File RelativePath="..\config.f90"/>
 		<File RelativePath="..\constants.f90"/>
-		<File RelativePath="..\cosmorec.f90"/>
 		<File RelativePath="..\DarkAge21cm.f90"/>
 		<File RelativePath="..\DarkEnergyFluid.f90"/>
 		<File RelativePath="..\DarkEnergyInterface.f90"/>
 		<File RelativePath="..\DarkEnergyPPF.f90"/>
 		<File RelativePath="..\equations.f90"/>
 		<File RelativePath="..\halofit.f90"/>
-		<File RelativePath="..\hyrec.f90"/>
 		<File RelativePath="..\InitialPower.f90"/>
 		<File RelativePath="..\lensing.f90"/>
 		<File RelativePath="..\massive_neutrinos.f90"/>

fortran/camb.f90

@@ -216,28 +216,6 @@ logical function CAMB_ReadParamFile(P, InputFile, InpLen, ErrMsg, ErrLen)
 
     end function CAMB_ReadParamFile
 
-    function UpperCase(str) Result (string)
-    !https://stackoverflow.com/questions/10759375/how-can-i-write-a-to-upper-or-to-lower-function-in-f90
-    !   ==============================
-    !   Changes a string to upper case
-    !   ==============================
-    Character(*), Intent(In) :: str
-    Character(LEN(str))      :: string
-
-    Integer :: ic, i
-
-    Character(26), Parameter :: cap = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
-    Character(26), Parameter :: low = 'abcdefghijklmnopqrstuvwxyz'
-
-    !   Capitalize each letter if it is lowecase
-    string = str
-    do i = 1, LEN_TRIM(str)
-        ic = INDEX(low, str(i:i))
-        if (ic > 0) string(i:i) = cap(ic:ic)
-    end do
-
-    End Function UpperCase
-
     logical function CAMB_ReadParams(P, Ini, ErrMsg)
     use NonLinear
     use DarkEnergyFluid

fortran/cmbmain.f90

@@ -152,7 +152,7 @@ subroutine cmbmain
 
     if (DebugMsgs .and. Feedbacklevel > 0) then
         call Timer%WriteTime('Timing for InitVars')
-        write (*,*) 'r = ',real(State%curvature_radius),' scale = ',real(State%scale), 'age = ', real(State%tau0)
+        if (.not. State%flat) call WriteFormat('r = %f, scale = %f',State%curvature_radius, State%scale)
     end if
 
     if (.not. State%OnlyTransfer .or. CP%NonLinear==NonLinear_Lens .or. CP%NonLinear==NonLinear_both) &
@@ -169,7 +169,7 @@ subroutine cmbmain
     !$ if (ThreadNum /=0) call OMP_SET_NUM_THREADS(ThreadNum)
 
     if (CP%WantCls) then
-        if (DebugMsgs .and. Feedbacklevel > 0) write(*,*) 'Set ',Evolve_q%npoints,' source k values'
+        if (DebugMsgs .and. Feedbacklevel > 0) call WriteFormat('Set %d source k values',Evolve_q%npoints)
 
         call GetSourceMem
 
@@ -230,7 +230,7 @@ subroutine cmbmain
 
             call SetkValuesForInt
 
-            if (DebugMsgs .and. Feedbacklevel > 0) write(*,*) 'Set ',ThisCT%q%npoints,' integration k values'
+            if (DebugMsgs .and. Feedbacklevel > 0) call WriteFormat('Set %d integration k values',ThisCT%q%npoints)
 
             !Begin k-loop and integrate Sources*Bessels over time
             !$OMP PARALLEL DO DEFAULT(SHARED), SCHEDULE(STATIC,4)
@@ -316,7 +316,7 @@ subroutine CalcLimberScalCls(CTrans)
     if (CP%SourceTerms%limber_phi_lmin>0) then
         winmin = 0
     else
-        winmin=1
+        winmin = 1
     end if
 
     do i =winmin, State%num_redshiftwindows
@@ -407,7 +407,8 @@ subroutine GetLimberTransfers
                 ThisCT%limber_l_min(s_ix) =  ell
                 ell_needed = ThisCT%ls%l(ell)
                 max_bessels_l_index = max(max_bessels_l_index,ThisCT%limber_l_min(s_ix)-1)
-                if (FeedbackLevel > 1) write (*,*) i,'Limber switch', ell_needed
+                if (FeedbackLevel > 1 .or. DebugMsgs) &
+                    call WriteFormat('Limber switch %d: %d', i,  ell_needed)
                 exit
             end if
         end do
@@ -476,7 +477,7 @@ subroutine GetLimberTransfers
                 end do
             end do
         else
-            max_bessels_l_index  = ThisCT%ls%nl
+            max_bessels_l_index = ThisCT%ls%nl
         end if
     end do
 
@@ -785,8 +786,7 @@ subroutine InitVars(state)
         end associate
     end do
 
-    !     Calculating the times for the outputs of the transfer functions.
-    !
+    ! Calculating the times for the outputs of the transfer functions.
     if (CP%WantTransfer .or. .not. allocated(State%Transfer_times)) then
         if (allocated(State%Transfer_times)) deallocate(State%Transfer_times)
         !allocate even if no transfer to prevent debugging access errors
@@ -1067,7 +1067,7 @@ subroutine TransferOut
 
 
     if (DebugMsgs .and. Feedbacklevel > 0) &
-        write(*,*) State%MT%num_q_trans-Evolve_q%npoints, 'transfer k values'
+        call WriteFormat('Transfer k values: %f',State%MT%num_q_trans-Evolve_q%npoints)
 
     !     loop over wavenumbers.
     !$OMP PARALLEL DO DEFAULT(SHARED),SCHEDULE(DYNAMIC), PRIVATE(EV, tau, q_ix)
@@ -1492,7 +1492,7 @@ subroutine DoFlatIntegration(IV, llmax)
             !Do integral if any useful contribution to the CMB, or large scale effects
 
             if (DoInt) then
-                if (CP%CUstomSources%num_custom_sources==0 .and. State%num_redshiftwindows==0) then
+                if (CP%CustomSources%num_custom_sources==0 .and. State%num_redshiftwindows==0) then
                     do n= State%TimeSteps%IndexOf(tmin),min(IV%SourceSteps,State%TimeSteps%IndexOf(tmax))
                         !Full Bessel integration
                         a2=aa(n)
@@ -1533,7 +1533,6 @@ subroutine DoFlatIntegration(IV, llmax)
                                 end do
                             end if
                         end do
-
                     else
                         do n= State%TimeSteps%IndexOf(tmin),min(IV%SourceSteps,State%TimeSteps%IndexOf(tmax))
                             !Full Bessel integration
@@ -2205,7 +2204,6 @@ subroutine CalcScalCls(CTrans)
     do j=1,CTrans%ls%nl
         !Integrate dk/k Delta_l_q**2 * Power(k)
         ell = real(CTrans%ls%l(j),dl)
-
         if (j<= CTrans%max_index_nonlimber) then
             do q_ix = 1, CTrans%q%npoints
                 if (.not.(State%closed.and.nint(CTrans%q%points(q_ix)*State%curvature_radius)<=CTrans%ls%l(j))) then
@@ -2334,8 +2332,6 @@ subroutine CalcScalCls(CTrans)
     !$OMP END PARALLEL DO
 #endif
 
-    deallocate(ks,pows,dlnks)
-
     end subroutine CalcScalCls
 
     subroutine CalcScalCls2(CTrans)

fortran/equations.f90

@@ -2896,7 +2896,6 @@ subroutine derivsv(EV,n,tau,yv,yvprime)
     !  sigma = 2*rhoq/k**2
     !for non-large k this expression for sigma is unstable at early times
     !so propagate sigma equation separately (near total cancellation in rhoq)
-    ! print *,yv(2),2*rhoq/k**2
 
     if (finished_tightcoupling) then
         !  Use explicit equations:

fortran/lensing.f90

@@ -738,7 +738,7 @@ subroutine BadHarmonic(State)
         if (RR > 1e-5) then
             write (*,*) 'You need to normalize realistically to use lensing.'
             write (*,*) 'see http://cosmocoffee.info/viewtopic.php?t=94'
-            call MpiStop()
+            call MpiStop('Lensing error')
         end if
         if (maxl > lmax_donelnfa) then
             !Get ln factorials

fortran/recfast.f90

@@ -978,8 +978,6 @@ subroutine ION(this,Ndim,z,Y,f)
             * (Tmat-Trad) / (Hz*(1._dl+z)) + 2._dl*Tmat/(1._dl+z)
     end if
 
-    ! print *, z, f(3)*(1+z)/Tmat
-
     if (evolve_Ts) then
 
         !       follow the matter temperature once it has a chance of diverging

fortran/results.f90

@@ -834,7 +834,6 @@ function CAMBdata_CosmomcTheta(this)
     rs = Integrate_Romberg(this,dsound_da_approx,1d-8,astar,atol)
     DA = this%AngularDiameterDistance(zstar)/astar
     CAMBdata_CosmomcTheta = rs/DA
-    !       print *,'z* = ',zstar, 'r_s = ',rs, 'DA = ',DA, rs/DA
 
     end function CAMBdata_CosmomcTheta
 
@@ -1767,7 +1766,8 @@ subroutine Thermo_Init(this, State,taumin)
                             RedWin%tau_start = tau01
                         else if (RedWin%tau_start /=0 .and. RedWin%tau_end==State%tau0 .and. winamp < 1e-8) then
                             RedWin%tau_end = min(State%tau0,tau + dtau)
-                            if (DebugMsgs) print *,'time window = ', RedWin%tau_start, RedWin%tau_end
+                            if (DebugMsgs) call WriteFormat('Window %d: tau1 = %f, tau2 = %f',&
+                                RW_i, RedWin%tau_start, RedWin%tau_end)
                         end if
                     else
                         if (RedWin%kind == window_lensing .or.  RedWin%kind == window_counts &
@@ -2143,7 +2143,7 @@ subroutine SetTimeSteps(this,State,TimeSteps)
     integer nri0, nstep
     !Sources
     integer ix,i,nwindow, L_limb
-    real(dl) keff, win_end, TimeSampleBoost
+    real(dl) keff, win_end, TimeSampleBoost, delta
 
     TimeSampleBoost = State%CP%Accuracy%AccuracyBoost*State%CP%Accuracy%TimeStepBoost
     call TimeSteps%Init()
@@ -2172,7 +2172,6 @@ subroutine SetTimeSteps(this,State,TimeSteps)
             associate (Win => State%Redshift_W(ix))
 
                 Win%tau_start = max(Win%tau_start,state%taurst)
-                this%tau_start_redshiftwindows = min(Win%tau_start,this%tau_start_redshiftwindows)
 
                 if (Win%kind /= window_lensing) then
                     !Have to be careful to integrate dwinV as the window tails off
@@ -2192,7 +2191,15 @@ subroutine SetTimeSteps(this,State,TimeSteps)
 
                 !Keep sampling in x better than Nyquist
                 nwindow = max(nwindow, nint(TimeSampleBoost *(win_end- Win%tau_start)* keff/3))
-                if (Feedbacklevel > 1) write (*,*) ix, 'nwindow =', nwindow
+                if (Win%kind /= window_lensing .and. Win%sigma_tau < (win_end - Win%tau_start)/nwindow*8) then
+                    nwindow = nint(TimeSampleBoost*(win_end - Win%tau_start)/Win%sigma_tau*8)
+                end if
+                if (Feedbacklevel > 1 .or. DebugMsgs) call WriteFormat('nwindow %d: %d', ix, nwindow)
+                delta = (win_end-Win%tau_start)/nwindow
+                Win%tau_start = Win%tau_start - delta*3
+                win_end = min(State%tau0, win_end+delta*2)
+
+                this%tau_start_redshiftwindows = min(Win%tau_start,this%tau_start_redshiftwindows)
 
                 call TimeSteps%Add(Win%tau_start, win_end, nwindow)
                 !This should cover whole range where not tiny
@@ -2243,7 +2250,7 @@ subroutine SetTimeSteps(this,State,TimeSteps)
         end do
     end if
 
-    if (DebugMsgs .and. FeedbackLevel > 0) write(*,*) 'Set ',nstep, ' time steps'
+    if (DebugMsgs .and. FeedbackLevel > 0) call WriteFormat('Set %d time steps', nstep)
 
     end subroutine SetTimeSteps
 
@@ -2653,7 +2660,7 @@ function Win_Limber_ell(W,CP,lmax) result(ell_limb)
         if (W%kind==window_lensing) then
             ell_limb = max(CP%SourceTerms%limber_phi_lmin,nint(50*LimBoost))
         else
-            ell_limb = max(CP%SourceTerms%limber_phi_lmin, nint(LimBoost *6* W%chi0/W%sigma_tau))
+            ell_limb = max(CP%SourceTerms%limber_phi_lmin, nint(LimBoost*6*W%chi0/W%sigma_tau))
         end if
     else
         ell_limb = lmax
@@ -3836,7 +3843,6 @@ subroutine Transfer_Get21cmCls(MTrans, State,FileNames)
     do itf_PK=1, CP%Transfer%PK_num_redshifts
         itf = State%PK_redshifts_index(itf_PK)
         if (FileNames(itf) /= '') then
-            !print *, 'tau = ', MTrans%optical_depths(itf)
             chi = State%tau0-State%TimeOfz(CP%Transfer%PK_redshifts(itf_PK))
 
             points = MTrans%num_q_trans

fortran/subroutines.f90

@@ -920,9 +920,7 @@ subroutine dverk (EV,n, fcn, x, y, xend, tol, ind, c, nw, w)
     !
 
     write (*,*) 'Error in dverk, x =',x, 'xend=', xend
-    call MpiStop()
-    !
-    !  end abort action
+    call MpiStop('DVERK error')
     !
     end subroutine dverk