anchor.diff:          nz=size(domain%u,2)
   
   !       load low-res U data
  -        call io_read3d(filename,options%uvar,extra_data,curstep)
  +        call io_read(filename,options%uvar,extra_data,curstep)
           domain%u=sum(extra_data(:,:,:nz))/size(extra_data(:,:,:nz))
           deallocate(extra_data)
   
   !       load low-res V data
  -        call io_read3d(filename,options%vvar,extra_data,curstep)
  +        call io_read(filename,options%vvar,extra_data,curstep)
           domain%v=sum(extra_data(:,:,:nz))/size(extra_data(:,:,:nz))
           deallocate(extra_data)
   
  
positive.code:         endif
  
          ! Read the data in, should be relatively fast because we are reading a low resolution forcing file
          call io_read3d(filename,varname,inputdata,curstep)
  
          nx=size(inputdata,1)
          ny=size(inputdata,2)
  
negative.code:         real,dimension(:,:),allocatable :: inputdata
  
  !       Read the data in
          call io_read(filename,varname,inputdata,curstep)
  !       interpolate data onto the high resolution grid
          call geo_interp2d(highres,inputdata,geolut)
          deallocate(inputdata)
  