anchor.diff:  
                       call interpolate_variable(var_to_interpolate%dqdt_3d, input_data, forcing, this, &
                                       vert_interp=var_is_not_pressure, var_is_u=var_is_u, var_is_v=var_is_v, nsmooth=this%nsmooth)
  +
  +                    ! because pressure needs to be adjusted for grid points that fall below the forcing lowest level, we adjust it separately.
                       if (.not.var_is_not_pressure) then
  -                        nz = min(size(this%geo%z, 2), size(forcing%geo%z, 2))
  -                        call update_pressure(var_to_interpolate%dqdt_3d, forcing%geo%z(:,:nz,:), this%geo%z)
  +                        allocate(potential_temperature, mold=var_to_interpolate%dqdt_3d)
  +                        ! to improve the pressure adjustment, we need to get forcing potential temperature on the ICAR grid WITHOUT vertical interpolation
  +                        call interpolate_variable(potential_temperature, forcing_temperature, forcing, this, &
  +                                        vert_interp=.False., var_is_u=.False., var_is_v=.False., nsmooth=this%nsmooth)
  +
  +                        call adjust_pressure(var_to_interpolate%dqdt_3d, forcing%geo%z, this%geo%z, potential_temperature)
                       endif
   
                   else
                       call interpolate_variable(var_to_interpolate%data_3d, input_data, forcing, this, &
                                       vert_interp=var_is_not_pressure, var_is_u=var_is_u, var_is_v=var_is_v, nsmooth=this%nsmooth)
   
  +                    ! because pressure needs to be adjusted for grid points that fall below the forcing lowest level, we adjust it separately.
                       if (.not.var_is_not_pressure) then
  -                        nz = min(size(this%geo%z, 2), size(forcing%geo%z, 2))
  -                        call update_pressure(var_to_interpolate%data_3d, forcing%geo%z(:,:nz,:), this%geo%z)
  +                        allocate(potential_temperature, mold=var_to_interpolate%dqdt_3d)
  +                        ! to improve the pressure adjustment, we need to get forcing potential temperature on the ICAR grid WITHOUT vertical interpolation
  +                        call interpolate_variable(potential_temperature, forcing_temperature, forcing, this, &
  +                                        vert_interp=.False., var_is_u=.False., var_is_v=.False., nsmooth=this%nsmooth)
  +
  +                        call adjust_pressure(var_to_interpolate%data_3d, forcing%geo%z, this%geo%z, potential_temperature)
                       endif
                   endif
   
  
positive.code: 
      end subroutine
  
  
      !> -------------------------------
      !! Interpolate one variable by requesting the forcing data from the boundary data structure then
  
negative.code:         ! temporary to hold the variable to be interpolated to
          type(variable_t) :: var_to_interpolate
          ! temporary to hold the forcing variable to be interpolated from
          type(variable_t) :: input_data, forcing_temperature
          real, allocatable, dimension(:,:,:) :: potential_temperature
          ! number of layers has to be used when subsetting for update_pressure (for now)
          integer :: nz
          logical :: var_is_u, var_is_v
  
          update_only = .False.
          if (present(update)) update_only = update
  
          forcing_temperature = forcing%variables%get_var(this%potential_temperature%meta_data%forcing_var)
  
          ! make sure the dictionary is reset to point to the first variable
          call this%variables_to_force%reset_iterator()
  