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:     end subroutine diagnostic_update
  
  
      !>------------------------------------------------------------
      !!  Uses the dt from step() to convert the forcing increments to per/time step increments
      !!
      !!  Divides dXdt variables by the length of the input timestep so that it can be multiplied by the internal
      !!  timestep to get the each update value.
      !!
      !! @param bc        Boundary conditions structure containing dXdt variables
      !! @param dt        Length of time between input forcing data for which the dXdt updates were calculated
      !! @param options   Model options structure specifies which variables need to be updated
      !!
      !!------------------------------------------------------------
      ! 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
      !
      !     ny=size(bc%du_dt,3)
      !     nx=size(bc%du_dt,1)
      !
      !     bc%du_dt  = bc%du_dt  / dt
      !     bc%dv_dt  = bc%dv_dt  / dt
      !     bc%dp_dt  = bc%dp_dt  / dt
      !     bc%dth_dt = bc%dth_dt / dt
      !     bc%dqv_dt = bc%dqv_dt / dt
      !     bc%dqc_dt = bc%dqc_dt / dt
      !
      !     if (trim(options%rain_var)/="") then
      !         bc%drain_dt = bc%drain_dt / dt
      !     endif
      !
      !     ! these only get updated if we are using the fluxes derived from the forcing model
      !     if (options%physics%landsurface==kLSM_BASIC) then
      !         bc%dsh_dt   = bc%dsh_dt   / dt
      !         bc%dlh_dt   = bc%dlh_dt   / dt
      !     endif
      !     ! these only get updated if we are using the fluxes derived from the forcing model
      !     if (options%physics%boundarylayer==kPBL_BASIC) then
      !         bc%dpblh_dt = bc%dpblh_dt / dt
      !     endif
      !     ! these only get updated if we are using the fluxes derived from the forcing model
      !     if (options%physics%radiation==kRA_BASIC) then
      !         bc%dsw_dt   = bc%dsw_dt   / dt
      !         bc%dlw_dt   = bc%dlw_dt   / dt
      !     endif
      !     bc%dsst_dt   = bc%dsst_dt   / dt
      !
      ! end subroutine apply_dt
      !
      !>------------------------------------------------------------
      !!  Calculate the maximum stable time step given some CFL criteria
      !!
  
negative.code:             ! step model_time forward
              domain%model_time = domain%model_time + dt
  
          enddo
  
      end subroutine step
  