anchor.diff:  ! ----------------------------------------------------------------------------
       use data_structures
       use io_routines,            only : io_getdims, io_read3d, io_maxDims, io_read2d, io_variable_is_present, &
  -                                       io_write3di, io_write3d
  +                                       io_write
       use wind,                   only : update_winds,balance_uvw
       use linear_theory_winds,    only : linear_perturb
       use geo,                    only : geo_interp2d, geo_interp
  
positive.code:         ! for pressure update we need real temperature, not potential t to compute an exner function
          ! bc%next_domain%pii=(bc%next_domain%p/100000.0)**(Rd/cp)
          ! now update pressure using the high res T field
          print*, "pre-update pressure Pmax",maxval(bc%next_domain%p), "pmin", minval(bc%next_domain%p)
          call update_pressure(bc%next_domain%p,bc%lowres_z,domain%z)!,  & 
                              !  lowresT = bc%next_domain%th * bc%next_domain%pii, &
                              !  hiresT  = domain%th * domain%pii)
          print*, "post-update pressure Pmax",maxval(bc%next_domain%p), "pmin", minval(bc%next_domain%p)
          
          ! not necessary as long as we are interpolating above
          call read_var(bc%next_domain%th,      file_list(curfile), options%tvar,   &
  
negative.code:             call geo_interp(inputdata, &
                              extra_data, &
                              geolut,boundary_only)
              ! Then apply vertical interpolation on that grid
              if (apply_vertical_interpolation) then
                  call vinterp(highres, inputdata, &
  