anchor.diff:  
   
       !>------------------------------------------------------------
  +    !!  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,
  
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:             ! step model_time forward
              domain%model_time = domain%model_time + dt
  
          enddo
  
      end subroutine step
  