anchor.diff:                    terrain               => this%terrain%data_2d)
   
               i = this%grid%kms
  -            ! print*, shape(dz_mass)
  -            ! print*, shape(dz)
  -            ! print*, shape(dz_interface)
  -            ! print*, shape(z_interface)
  -            ! print*, shape(z)
  -            dz_mass(:,i,:) = dz(1) / 2
  -            z(:,i,:) = terrain + dz_mass(:,i,:)
  -            z_interface(:,i,:) = terrain
  -
  -            ! stop "need to think about dz indexing here"
  +            dz_mass(:,i,:)      = dz(i) / 2
  +            z(:,i,:)            = terrain + dz_mass(:,i,:)
  +            z_interface(:,i,:)  = terrain
  +
               do i = this%grid%kms+1, this%grid%kme
  -                dz_mass(:,i,:)     = (dz(i-this%grid%kms+1) + dz(i-this%grid%kms)) / 2
  -                dz_interface(:,i,:)= dz(i-this%grid%kms+1)
  +                dz_mass(:,i,:)     = (dz(i) + dz(i-1)) / 2
  +                dz_interface(:,i,:)= dz(i)
                   z(:,i,:)           = z(:,i-1,:)           + dz_mass(:,i,:)
                   z_interface(:,i,:) = z_interface(:,i-1,:) + dz_interface(:,i,:)
               enddo
  -            ! note interface elements need to be specified at an additional level
  -            ! i = this%grid%kme + 1
  -            ! z_interface(:,i,:) = z_interface(:,i-1,:) + dz_interface(:,i,:)
  -
  -            exner             = exner_function(pressure)
  -            temperature       = exner * potential_temperature
  +
           end associate
   
       end subroutine initialize_variables
  
positive.code: 
          ! List the variables that are required for any restart
          call options%restart_vars(                                                  &
                       [kVARS%pressure,                                               &
                        kVARS%potential_temperature,                                  &
                        kVARS%z,                                                      &
                        kVARS%terrain,                                                &
                        kVARS%latitude,               kVARS%longitude,                &
                        kVARS%u_latitude,             kVARS%u_longitude,              &
  
negative.code:       class(domain_t), intent(inout) :: this
        class(options_t),intent(inout) :: options
  
        call mp_simple_var_request(options)
  
        call this%var_request(options)
  
        call read_domain_shape(this, options)
  
        call create_variables(this, options)
  
  