anchor.diff:              model_time = model_time + dt
               domain%model_time = model_time
   
  -            if ((abs(model_time-next_output)<1e-1).or.(model_time>next_output)) then
  -                call write_domain(domain,options,nint((model_time-options%time_zero)/options%out_dt))
  +            if ((abs(model_time - next_output) < 1e-1).or.(model_time > next_output)) then
  +                call write_domain(domain, options, nint((next_output - options%time_zero)/options%out_dt))
  +                ! Note that over a long simulation, this can build up errors.
  +                ! Need to add some checks with respect to model initial time...
                   next_output = next_output + options%out_dt
               endif
   
  
positive.code:     subroutine apply_dt(bc, dt, options)
          implicit none
          type(bc_type),      intent(inout)   :: bc
          real,               intent(in)      :: dt
          type(options_type), intent(in)      :: options
          ! internal variables
          integer :: j, ny, nx
  
negative.code:         type(options_type), intent(in)      :: options
          real*8,             intent(inout)   :: model_time, next_output
  
          real :: time_percent
          real :: dt, CFL, cfl_reduction_factor
  
          CFL = 1.0
          cfl_reduction_factor = options%cfl_reduction_factor
          CFL = CFL * cfl_reduction_factor
  
          ! Make the boundary condition dXdt values into units of [X]/s
          call apply_dt(bc, bc%dt, options)
  
          ! ensure internal model consistency (should only need to be called here when the model starts...)
          ! e.g. for potential temperature and temperature
          call diagnostic_update(domain,options)
  
          ! now just loop over internal timesteps computing all physics in order (operator splitting...)
          do while (model_time < bc%next_domain%model_time)
  
              ! compute internal timestep dt to maintain stability
              ! courant condition for 3D advection. Note that w is normalized by dx/dz
  