anchor.diff:          type(Time_type),    intent(in)      :: end_time
           type(options_t),    intent(in)      :: options
   
  -        real :: time_percent, last_time
  -        double precision :: seconds
  -        type(time_delta_t) :: dt, progress_dt, time_step_size
  +        real :: last_print_time
  +        type(time_delta_t) :: dt, time_step_size
   
           call update_winds(domain, options)
   
  +        last_print_time = 0.0
           time_step_size = end_time - domain%model_time
  +
           ! Make the boundary condition dXdt values into units of [X]/s
           call domain%update_delta_fields(time_step_size)
   
   
           ! now just loop over internal timesteps computing all physics in order (operator splitting...)
  -        last_time = 0.0
           do while (domain%model_time < end_time)
   
               ! ensure internal model consistency
               call diagnostic_update(domain, options)
   
  -            ! compute internal timestep dt to maintain stability
  -            ! courant condition for 3D advection. Note that w is normalized by dx/dz
  -            ! pick the minimum dt from the begining or the end of the current timestep
  -            if (options%parameters%advect_density) then
  -                ! call dt%set(seconds=compute_dt(domain%dx, domain%ur, domain%vr, domain%wr, domain%rho, domain%dz_inter,&
  -                !                 options%parameters%cfl_reduction_factor, cfl_strictness=options%parameters%cfl_strictness,                   &
  -                !                 use_density=.True.))
  -
  -            else
  -                ! compute the dt to meet the CFL criteria specified given dx, u, v, w, dz
  -                call dt%set(seconds=compute_dt(domain%dx, domain%u%data_3d, domain%v%data_3d, domain%w%data_3d, domain%density%data_3d, domain%dz_interface%data_3d, &
  -                                options%parameters%cfl_reduction_factor, cfl_strictness=options%parameters%cfl_strictness,                 &
  -                                use_density=.False.))
  -            endif
  -
  -            ! set an upper bound on dt to keep microphysics and convection stable (?)
  -            seconds = dt%seconds()
  -            ! perform a reduction across all images to find the minimum time step required
  -            sync all ! not sure if this is required before the co_min call
  -            call co_min(seconds)
  -            ! store this back in the dt time_delta data structure
  -            call dt%set(seconds=min(seconds,120.0D0))
  -
  +            call update_dt(dt, options, domain, end_time)
   
               ! if an interactive run was requested than print status updates everytime at least 5% of the progress has been made
               if (options%parameters%interactive .and. (this_image()==1)) then
  -
  -                ! first compute the current time until reaching the end
  -                progress_dt  = (end_time - domain%model_time)
  -                ! convert that to a percentage of the total time required
  -                time_percent = 100 - progress_dt%seconds() / time_step_size%seconds()  * 100
  -                ! finally if it has been at least 5% of the time since the last time we printed output, print output
  -                if (time_percent > (last_time + 5.0)) then
  -                    last_time = last_time + 5.0
  -                    ! this used to just use the nice $ (or advance="NO") trick, but at least with some mpi implementations, it buffers this output until it crashes
  -                    write(*,"(A,f5.1,A,A)") char(13), max(0.0,time_percent)," %  dt=",trim(dt%as_string())
  -                endif
  -            endif
  -
  -            ! Make sure we don't over step the forcing period
  -            if ((domain%model_time + dt) > end_time) then
  -                dt = end_time - domain%model_time
  +                call print_progress(domain%model_time, end_time, time_step_size, dt, last_print_time)
               endif
   
               ! this if is to avoid round off errors causing an additional physics call that won't really do anything
  
positive.code:         type(domain_t),  intent(inout)   :: domain
          type(options_t), intent(in)      :: options
  
          domain%exner%data_3d = exner_function(domain%pressure%data_3d)
  
          domain%density%data_3d =  domain%pressure%data_3d / (Rd * domain%potential_temperature%data_3d * domain%exner%data_3d) ! kg/m^3
  
  
      !
      !     integer :: nx,ny,nz, y, z
      !
  
negative.code: 
  
      !>------------------------------------------------------------
      !!  Prints progress to the terminal if requested
      !!
      !! @param current_time  the current state of the model time
      !! @param end_time      the end of the current full time step (when step will return)
      !! @param time_step     length of full time step to be integrated by step
      !! @param dt            numerical timestep to print for information
      !!
      !!------------------------------------------------------------
      subroutine print_progress(current_time, end_time, time_step, dt, last_time)
          implicit none
          type(Time_type),    intent(in)    :: current_time,    end_time
          type(time_delta_t), intent(in)    :: time_step,       dt
          real,               intent(inout) :: last_time
  
          type(time_delta_t) :: progress_dt
          real :: time_percent
  
          ! first compute the current time until reaching the end
          progress_dt  = (end_time - current_time)
  
          ! convert that to a percentage of the total time required
          time_percent = 100 - progress_dt%seconds() / time_step%seconds()  * 100
  
          ! finally if it has been at least 5% of the time since the last time we printed output, print output
          if (time_percent > (last_time + 5.0)) then
              last_time = last_time + 5.0
              ! this used to just use the nice $ (or advance="NO") trick, but at least with some mpi implementations, it buffers this output until it crashes
              write(*,"(A,f5.1,A,A)") char(13), max(0.0, time_percent)," %  dt=",trim(dt%as_string())
          endif
  
      end subroutine print_progress
  
      !>------------------------------------------------------------
      !! Update the numerical timestep to use
      !!
      !! @param dt            numerical timestep to use
      !! @param options       set options for how to update the time step
      !! @param domain        the full domain structure (need winds to compute dt)
      !! @param end_time      the end of the current full time step (when step will return)
      !!
      !!------------------------------------------------------------
      subroutine update_dt(dt, options, domain, end_time)
          implicit none
          type(time_delta_t), intent(inout) :: dt
          type(options_t),    intent(in)    :: options
          type(domain_t),     intent(in)    :: domain
          type(Time_type),    intent(in)    :: end_time
  
          double precision :: seconds
  
          ! compute internal timestep dt to maintain stability
          ! courant condition for 3D advection. Note that w is normalized by dx/dz
  
          ! Note this needs to be performed when advect_density is enabled
          ! if (options%parameters%advect_density) then
              ! call dt%set(seconds=compute_dt(domain%dx, domain%ur, domain%vr, domain%wr, domain%rho, domain%dz_inter,&
              !                 options%parameters%cfl_reduction_factor, cfl_strictness=options%parameters%cfl_strictness,                   &
              !                 use_density=.True.))
  
          ! else
          ! compute the dt to meet the CFL criteria specified given dx, u, v, w, dz
          associate(dx         => domain%dx,                              &
                    u          => domain%u%data_3d,                       &
                    v          => domain%v%data_3d,                       &
                    w          => domain%w%data_3d,                       &
                    density    => domain%density%data_3d,                 &
                    dz         => domain%dz_interface%data_3d,            &
                    reduction  => options%parameters%cfl_reduction_factor,&
                    strictness => options%parameters%cfl_strictness       &
              )
  
              seconds = compute_dt(dx, u, v, w, density, dz, reduction, &
                                   cfl_strictness=strictness, use_density=.false.)
  
          end associate
          ! endif
  
          ! perform a reduction across all images to find the minimum time step required
          call co_min(seconds)
  
          ! set an upper bound on dt to keep microphysics and convection stable (?)
          ! store this back in the dt time_delta data structure
          call dt%set(seconds=min(seconds,120.0D0))
  
          ! Make sure we don't over step the forcing period
          if ((domain%model_time + dt) > end_time) then
              dt = end_time - domain%model_time
          endif
  
      end subroutine update_dt
  
  
      !>------------------------------------------------------------
      !!  Step forward one IO time step.
      !!
      !!  Calculated the internal model time step to satisfy the CFL criteria,
  