%% copied from /home/sergio/HITRAN2UMBCLBL/make_refprof_usstd2010_lbl.m

clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iV = 1;  %% July 1, 2010
iV = 2;  %% Aug 06, 2010 Scott updated profs for gas 41/80
iV = 3;  %% Aug 16, 2010 Scott updated profs for gas 41/80
if iV == 1
  load /asl/packages/klayersV205/Data/refprof_usstd2010.mat
elseif iV == 2
  load /asl/packages/klayersV205/Data/refprof_usstd6Aug2010.mat
elseif iV == 3
  load /asl/packages/klayersV205/Data/refprof_usstd16Aug2010.mat
  end

% p and pp must be in atmospheres;   see /SPECTRA/DOC/lbl.tex
% T        must be in Kelvin
% q        must be in kilomoles/cm2

for ii = 1 : 100
  num   = prof.plevs(ii)-prof.plevs(ii+1);
  denom = log(prof.plevs(ii)/prof.plevs(ii+1));
  mpres(ii) = (num/denom)/1013.25;   %% change from mb to atm
  mtemp(ii)  = prof.ptemp(ii);
  thickness(ii) = abs(prof.palts(ii) - prof.palts(ii+1));   %% in meters
  qall(ii) = mpres(ii)*101325 * thickness(ii)/(8.314674269981136 * mtemp(ii)); 
  %% change from moles/m2 to moles/cm2 to kilomoles/cm2
  qall(ii) = qall(ii)/1e4/1e3;
  end
plev = prof.plevs;

glistall = [[1:42] [51:81]];
for ii = 1 : length(glistall)
  gas = glistall(ii);
  glist(ii) = gas;
  str = ['thegas = prof.gas_' num2str(gas) ';']; eval(str);
  thegas = thegas(1:100);             %% in molecules/cm2
  gamnt(:,ii) = thegas/6.022045e23/1000; %% in kilomoles/cm2
  gpart(:,ii) = gamnt(:,ii)./qall' .* mpres';
  end

glist = glist';
mpres = mpres'; mpres = flipud(mpres);
mtemp = mtemp'; mtemp = flipud(mtemp);
refpro.glist = glist;
refpro.mpres = mpres;
refpro.mtemp = mtemp;
refpro.gamnt = flipud(gamnt);
refpro.gpart = flipud(gpart);
refpro.plev  = flipud(plev);

for ii = 1 : length(glistall)
  if refpro.glist(ii) ~= glistall(ii)
    [refpro.glist(ii)  glistall(ii)]
    error('gasids do not match!')
    end
  outfile = ['refgas' num2str(glistall(ii))];
  fid = fopen(outfile,'w');
  if iV == 1
    string = '!Ref Profile, 07/01/2010 /asl/packages/klayersV205/Data/refprof_usstd2010.mat';
  elseif iV == 2
    string = '!Ref Profile, 08/10/2010 /asl/packages/klayersV205/Data/refprof_usstd6Aug2010.mat';
  elseif iV == 3
    string = '!Ref Profile, 08/19/2010 /asl/packages/klayersV205/Data/refprof_usstd16Aug2010.mat';
    end
  fprintf(fid,'%s\n',string);
  string = ['! GasID ' num2str(glistall(ii))];
  fprintf(fid,'%s\n',string);
  string = '! lay    pressure    part pres    temp      amount';
  fprintf(fid,'%s\n',string);
  string = '!----  -----------  -----------  -------  -----------';
  fprintf(fid,'%s\n',string);
  data = [glistall(ii)*ones(100,1) ...
        refpro.mpres refpro.gpart(:,ii) refpro.mtemp refpro.gamnt(:,ii)];
  fprintf(fid,'%2i  %10.4e  %10.4e  %7.3f  %10.4e \n',data');
  fclose(fid);
  end

toucher = ['!touch output_refprof_usstdJuly2010.m'];
eval(toucher)