anchor.diff:          type(options_type), intent(in)      :: options
           real*8,             intent(inout)   :: model_time, next_output
   
  -        real*8  :: end_time
  -        real    :: dt, CFL, cfl_reduction_factor
  +        real :: time_percent
  +        real :: dt, CFL, cfl_reduction_factor
   
           CFL = 1.0
           cfl_reduction_factor = options%cfl_reduction_factor
           CFL = CFL * cfl_reduction_factor
   
  -        ! calculate the number of timesteps
  -        end_time = model_time + options%in_dt
  -
           ! Make the boundary condition dXdt values into units of [X]/s
  -        call apply_dt(bc, options%in_dt, options)
  +        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 < end_time)
  +        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
  
positive.code:             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))
                  next_output = next_output + options%out_dt
              endif
  
  
negative.code:     subroutine apply_dt(bc, dt, options)
          implicit none
          type(bc_type),      intent(inout)   :: bc
          double precision,   intent(in)      :: dt
          type(options_type), intent(in)      :: options
          ! internal variables
          integer :: j, ny, nx
  