anchor.diff:  
           call setup_geo(this%geo,   this%latitude%data_2d,   this%longitude%data_2d,   this%z%data_3d)
   
  -    end subroutine initialize_variables
  +    end subroutine initialize_core_variables
  +
  +    !> -------------------------------
  +    !! Initialize various internal variables that need forcing data first, e.g. temperature, pressure on interface, exner, ...
  +    !!
  +    !! -------------------------------
  +    subroutine initialize_internal_variables(this, options)
  +        implicit none
  +        class(domain_t), intent(inout)  :: this
  +        class(options_t),intent(in)     :: options
  +
  +        integer :: i
  +
  +        associate(pressure              => this%pressure%data_3d,               &
  +                  exner                 => this%exner%data_3d,                  &
  +                  pressure_interface    => this%pressure_interface%data_3d,     &
  +                  psfc                  => this%surface_pressure%data_2d,       &
  +                  temperature           => this%temperature%data_3d,            &
  +                  potential_temperature => this%potential_temperature%data_3d )
  +
  +                  exner = exner_function(pressure)
  +
  +                  if (associated(this%pressure_interface%data_3d)) then
  +                      ! this isn't exactly correct, should be distance weighted...
  +                      ! weight one = (dz2) / (dz1+dz2)
  +                      ! weight two = (dz1) / (dz1+dz2)
  +                      pressure_interface(:,1,:) = ( pressure(:,1,:) * 2 - pressure(:,2,:) )
  +                      do i = 2, size(pressure_interface, 2)
  +                          pressure_interface(:,i,:) = ( pressure(:,i-1,:) + pressure(:,i,:) ) / 2
  +                      enddo
  +
  +                      if (associated(this%surface_pressure%data_2d)) then
  +                          psfc = pressure_interface(:,1,:)
  +                      endif
  +                  endif
  +
  +                  if (associated(this%temperature%data_3d)) then
  +                      temperature = potential_temperature * exner
  +                  endif
  +
  +        end associate
  +
  +    end subroutine initialize_internal_variables
   
       !> -------------------------------
       !! Populare the metadata structure in the domain for later output
  
positive.code: 
          call create_variables(this, options)
  
          call initialize_variables(this, options)
  
          call setup_meta_data(this, options)
  
  
negative.code:     !! Initialize various domain variables, mostly z, dz, etc.
      !!
      !! -------------------------------
      subroutine initialize_core_variables(this, options)
          implicit none
          class(domain_t), intent(inout)  :: this
          class(options_t),intent(in)     :: options
  