anchor.diff:  
       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
  
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: 
                      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
                          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
                          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
  
  