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:         ! 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
          ! 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
  
          ! make sure the dictionary is reset to point to the first variable
          call this%variables_to_force%reset_iterator()
  
negative.code: 
      end subroutine
  
      !> -------------------------------
      !! Adjust a 3d pressure field from the forcing data to the ICAR model grid
      !!
      !! Because the GCM grid can be very different from the ICAR grid, we first roughly match up
      !! the GCM level that is closest to the ICAR level. This has to be done grid cell by gridcell.
      !! This still is not ideal, in that it has already subset the GCM levels to the same number as are in ICAR
      !! If the GCM has a LOT of fine layers ICAR will not be getting layers higher up in the atmosphere.
      !! It would be nice to first use vinterp to get as close as we can, then update pressure only for grid cells below.
      !! Uses update_pressure to make a final adjustment (including below the lowest model level).
      !!
      !! -------------------------------
      subroutine adjust_pressure(pressure, input_z, output_z, potential_temperature)
          implicit none
          real, intent(inout), dimension(:,:,:) :: pressure !> Pressure on the forcing model levels [Pa]
          real, intent(in), dimension(:,:,:) :: input_z, output_z !> z on the forcing and ICAR model levels [m]
          real, intent(in), dimension(:,:,:) :: potential_temperature !> potential temperature of the forcing data [K]
  
          ! store a temporary copy of P and Z from the forcing data after selecting the closest GCM level to the ICAR data
          real, allocatable, dimension(:,:,:) :: temp_z, temp_p, temp_t
          ! loop counter variables
          integer :: k, nz, in_z_idx
          integer :: i,j, nx, ny
  
          allocate(temp_z, temp_p, temp_t, mold=pressure)
  
          nx = size(pressure, 1)
          nz = size(pressure, 2)
          ny = size(pressure, 3)
  
          do j = 1, ny
              do i = 1, nx
                  ! keep track of the nearest z level from the forcing data
                  in_z_idx = 1
                  do k = 1, nz
                      ! if the ICAR z level is more than half way to the next forcing z level, then increment the GCM z
                      findz: do while (output_z(i,k,j) > ((input_z(i,in_z_idx,j) + input_z(i,min(nz,in_z_idx+1),j)) / 2))
                          in_z_idx = min(nz, in_z_idx + 1)
  
                          if (in_z_idx == nz) then
                              exit findz
                          endif
                      end do findz
                      ! make a new copy of the pressure and z data from the closest GCM model level
                      temp_z(i,k,j) = input_z(i,in_z_idx,j)
                      temp_p(i,k,j) = pressure(i,in_z_idx,j)
                      temp_t(i,k,j) = exner_function(pressure(i,in_z_idx,j)) * potential_temperature(i,in_z_idx,j)
                  end do
              enddo
          enddo
  
          ! put the updated pressure data into the pressure variable prior to adjustments
          pressure = temp_p
  
          ! update pressure for the change in height between the closest GCM model level and each ICAR level.
          call update_pressure(pressure, temp_z, output_z, temp_t)
  
          deallocate(temp_p, temp_z)
      end subroutine
  
      !> -------------------------------
      !! Interpolate one variable by requesting the forcing data from the boundary data structure then
  