function jac_rtp(fin, fout, iprof, dT, dA, glist);

% function jac_rtp(fin, fout, iprof, dT, dA, glist);
%
% Create an RTP file for doing finite difference Jacobians.
%
% Input:
%    fin = (string) name of RTP input file
%    fout = (string) name of RTP output file
%    iprof = (1 x 1) profile number (index in input file)
%    dT = (1 x 1) temperature perturbation, in Kelvin (typically 1)
%    dA = (1 x 1) gas amount perturbation, in percent (typically 10)
%    glist = (ngas x 1) gas ids for amount perturbation
%
% Output: none (see output file)
%

% Created 12 September 2001, Scott Hannon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ngjac=length(glist);

addpath /asl/matlab/h4tools

% Read input file
[headin, hattr, profin, pattr]=rtpread(fin);


%%%%%%%%%%%%%%
% Check header
%%%%%%%%%%%%%%
if (headin.ptype < 1)
   disp('ERROR: input file is not a "layers" profile');
   return
end
%
if (ngjac > 0)
   [j1, j2, gindex]=intersect(glist, headin.glist);
   if (length(gindex) ~= ngjac)
      disp('Error: input file missing some gases for jacobians');
      disp('Selected gases:')
      glist
      disp('Found gases:')
      headin.glist(gindex)
      return
   end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create output file "head" structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
head.ptype=headin.ptype;
head.pfields=1;
head.ngas=headin.ngas;
head.glist=headin.glist;
head.gunit=headin.gunit;
if (isfield(headin,'nchan'))
   head.nchan=headin.nchan;
end
if (isfield(headin,'ichan'))
   head.ichan=headin.ichan;
end
head.mwnchan=0;
%
clear headin


% Check profile number
if (iprof > length(profin.nlevs))
   disp('ERROR: bad profile number')
   return
else
   nlevs=profin.nlevs(iprof);
   nlay=nlevs - 1;
end


% Determine number of profile entries for output file
% nout=unperturbed + nlay*(temperature + each selected gas)
nout=1 + nlay*(1 + ngjac);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create output "prof" structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Location
if (isfield(profin,'plat'))
   prof.plat=profin.plat(iprof)*ones(1,nout);
end
if (isfield(profin,'plon'))
   prof.plon=profin.plon(iprof)*ones(1,nout);
end
if (isfield(profin,'ptime'))
   prof.ptime=profin.ptime(iprof)*ones(1,nout);
end
if (isfield(profin,'landfrac'))
   prof.landfrac=profin.landfrac(iprof)*ones(1,nout);
end
%
% Surface
if (isfield(profin,'stemp'))
   prof.stemp=profin.stemp(iprof)*ones(1,nout);
end
if (isfield(profin,'spres'))
   prof.spres=profin.spres(iprof)*ones(1,nout);
end
if (isfield(profin,'salti'))
   prof.salti=profin.salti(iprof)*ones(1,nout);
end
%
% Emissivity & reflectivity
if (isfield(profin,'nrho'))
   nrho=profin.nrho(iprof);
   prof.nrho=nrho*ones(1,nout);
   prof.rfreq=profin.rfreq(1:nrho,iprof)*ones(1,nout);
   prof.rho=profin.rho(1:nrho,iprof)*ones(1,nout);
end
if (isfield(profin,'nemis'))
   nemis=profin.nemis(iprof);
   prof.nemis=nemis*ones(1,nout);
   prof.efreq=profin.efreq(1:nemis,iprof)*ones(1,nout);
   prof.emis=profin.emis(1:nemis,iprof)*ones(1,nout);
end
%
% Angles
if (isfield(profin,'scanang'))
   prof.scanang=profin.scanang(iprof)*ones(1,nout);
end
if (isfield(profin,'satzen'))
   prof.satzen=profin.satzen(iprof)*ones(1,nout);
end
if (isfield(profin,'satazi'))
   prof.satazi=profin.satazi(iprof)*ones(1,nout);
end
if (isfield(profin,'solzen'))
   prof.solzen=profin.solzen(iprof)*ones(1,nout);
end
if (isfield(profin,'solazi'))
   prof.solazi=profin.solazi(iprof)*ones(1,nout);
end
%
% Profiles
prof.nlevs=nlevs*ones(1,nout);
if (isfield(profin,'plevs'))
   prof.plevs=profin.plevs(1:nlevs,iprof)*ones(1,nout);
end
if (isfield(profin,'plays'))
   prof.plays=profin.plays(1:nlevs,iprof)*ones(1,nout);
end
if (isfield(profin,'palts'))
   prof.palts=profin.palts(1:nlevs,iprof)*ones(1,nout);
end
prof.ptemp=profin.ptemp(1:nlevs,iprof)*ones(1,nout);
for ig=1:head.ngas
   gstr=int2str(head.glist(ig));
   eval(['prof.gas_' gstr '=profin.gas_' gstr '(1:nlevs,iprof)*ones(1,nout);']);
end
if (isfield(profin,'co2ppm'))
   prof.co2ppm=profin.co2ppm(iprof)*ones(1,nout);
end
%
% No clouds, no microwave, no observations
%
clear profin


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjust prof for the perturbations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Temperature
ip=2;
for il=1:nlay
   prof.ptemp(il,ip)=prof.ptemp(il,ip) + dT;
   ip=ip + 1;
end
%
% Amount
x=1 + (dA/100);
for ig=1:ngjac
   gstr=int2str( head.glist(gindex(ig)) );
   for il=1:nlay
      eval(['prof.gas_' gstr '(il,ip)=prof.gas_' gstr '(il,ip)*x;']);
      ip=ip + 1;
   end
end


%%%%%%%%%%%%%%%%%%%
% Write output file
%%%%%%%%%%%%%%%%%%%
rtpwrite(fout, head, hattr, prof, pattr)


%%% end of file %%%
