anchor.diff:          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
  -
  -
  +        associate(exner                 => domain%exner%data_3d,    &
  +                  pressure              => domain%pressure%data_3d, &
  +                  density               => domain%density%data_3d,  &
  +                  potential_temperature => domain%potential_temperature%data_3d)
  +
  +        exner = exner_function(pressure)
  +
  +        density =  pressure / &
  +                    (Rd * potential_temperature * exner) ! kg/m^3
  +
  +        end associate
  +
  +    ! NOTE: all code below is not implemented in ICAR 2.0 yet
  +    ! it is left as a reminder of what needs to be done, and example when the time comes
       !
       !     integer :: nx,ny,nz, y, z
       !
  
positive.code:             ! step model_time forward
              domain%model_time = domain%model_time + dt
  
              progress_dt = domain%model_time - end_time ! note, this was next_output...
  
              ! if ((abs(progress_dt%seconds()) < 1e-1).or.(domain%model_time > next_output)) then
              !     write(*,*) " "
              !     call write_domain(domain, options, next_output)
              !     ! 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%output_dt
              ! endif
  
          enddo
  
      end subroutine step
  
negative.code:         type(Time_type),    intent(in)      :: end_time
          type(options_t),    intent(in)      :: options
  
          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...)
          do while (domain%model_time < end_time)
  
              ! ensure internal model consistency
              call diagnostic_update(domain, options)
  
              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
                  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
  