function [y,J] = arts_oem( Q, R, x, iter, do_j )


%- Set default values
%
do_j_DEFAULT = 0;
%
set_defaults;

Q0=Q;

% Checks done *arts_oem_init* are not repeated here


%--- Map x to variables -------------------------------------------------------
if iter > 1

  %- Some other variables
  %
  i_asj = find( [ Q.ABS_SPECIES_JAC.DO ] );
  nq    = length( R.jq );

  %- Loop retrieval quantities
  %
  for i = 1 : nq

    ind = R.ji{i}{1} : R.ji{i}{2};

    switch R.jq{i}.maintag

    case 'Gas species'   %---------------------------------------------------
      %
      ig  = i_asj(i);    % Gas species index
      D   = [];
      %
      if strcmp( R.jq{i}.mode, 'rel' )
        [D.P_GRID,D.LAT_GRID,D.LON_GRID,X] = x2field( Q.ATMOSPHERE_DIM, ...
            Q.ABS_SPECIES_JAC(ig).P_GRID, Q.ABS_SPECIES_JAC(ig).LAT_GRID, ...
            Q.ABS_SPECIES_JAC(ig).LON_GRID, x(ind) .* R.vmr(ind), ...
                                            Q.P_GRID, Q.LAT_GRID, Q.LON_GRID );
      elseif strcmp( R.jq{i}.mode, 'vmr' )
        [D.P_GRID,D.LAT_GRID,D.LON_GRID,X] = x2field( Q.ATMOSPHERE_DIM, ...
            Q.ABS_SPECIES_JAC(ig).P_GRID, Q.ABS_SPECIES_JAC(ig).LAT_GRID, ...
            Q.ABS_SPECIES_JAC(ig).LON_GRID, x(ind), ...
                                            Q.P_GRID, Q.LAT_GRID, Q.LON_GRID );
      else
        [D.P_GRID,D.LAT_GRID,D.LON_GRID,X] = x2field( Q.ATMOSPHERE_DIM, ...
           Q.ABS_SPECIES_JAC(ig).P_GRID, Q.ABS_SPECIES_JAC(ig).LAT_GRID, ...
           Q.ABS_SPECIES_JAC(ig).LON_GRID, nd2vmr(x(ind),R.p(ind),R.t(ind)),...
                                            Q.P_GRID, Q.LAT_GRID, Q.LON_GRID );
      end
      %
      Q.VMR_FIELD(ig,:,:,:) = arts_atminterp( Q.ATMOSPHERE_DIM, D, X, Q );
      %
    end
  end
end
%------------------------------------------------------------------------------



%--- Run ARTS -----------------------------------------------------------------

if nargout == 1  |  ~do_j
  y = arts_y( Q );
else
  [y,dy,J] = arts_y( Q );
end
%
%------------------------------------------------------------------------------


return





%--- Internal sub-functions ---------------------------------------------------

function [rgrid1,rgrid2,rgrid3,x] = ...
                 x2field( dim, rgrid1, rgrid2, rgrid3, x, grid1, grid2, grid3 )
  %
  if dim == 1
    %
    [rgrid1,f1,l1] = expand_grid( rgrid1, grid1 );
    %
    if ~isempty(f1)  |  ~isempty(l1)
      x = [ x(f1); vec2col(x); x(l1) ];   % x as row vector accepted !
    end 
    %
  elseif dim == 2
    %
    x = reshape( x, length(rgrid1), length(rgrid2) );
    %
    [rgrid1,f1,l1] = expand_grid( rgrid1, grid1 );
    [rgrid2,f2,l2] = expand_grid( rgrid2, grid2 );
    %
    if ~isempty(f1)  |  ~isempty(l1)  |  ~isempty(f2)  |  ~isempty(l2)
      x = [ x(f1,f2) x(f1,:) x(f1,l2); 
            x(:,f2)  x       x(:,l2); 
            x(l1,f2) x(l1,:) x(l1,l2) ];
    end 
  end
return


%------------------------------------------------------------------------------

function [rgrid,f,l] = expand_grid( rgrid, grid )
  %
  % *rgrid* is expanded to cover range of *grid*. The output arguments
  % *f* and *l* can be used to expand the data to be interpolated following
  % the usage at the end of this function.
  %
  n1 = length( grid );
  n2 = length( rgrid );
  %
  % Increasing vector
  if grid(2) > grid(1)
    if grid(1) < rgrid(1)
      f = 1;
    else
      f = [];
    end 
    if grid(n1) > rgrid(n2)
      ll = n1; l = n2;
    else
      ll = []; l = [];
    end 
  % Decreasing vector
  else
    if grid(1) > rgrid(1)
      f = 1;
    else
      f = [];
    end 
    if grid(n1) < rgrid(n2)
      ll = n1; l = n2;
    else
      ll = []; l = [];
    end 
  end
  %
  if ~isempty(f)  |  ~isempty(l)
    rgrid = [ grid(f); vec2col(rgrid); grid(ll) ];
  end
  %
return