%%%%AIRS levels
A = -1.5508E-4;
B = -5.5937E-2;
C =  7.4516;
j = 1 : 101;
Plev = exp( (7/2) * log( A*j.^2 + B*j + C ) );
plot(j,Plev)

husual = p2h(Plev)/1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
caIn = ['trp','mls','mlw','sas','saw','std'];

%%plot the (T,p) profs
for iI = 1 : 6
  iInd = (1:3) + (iI-1)*3;
  caX = caIn(iInd);
  prof = load(caX);
  semilogy(prof(:,4),prof(:,3)); hold on
  end
%check that US std +/- 50 covers everything
semilogy(prof(:,4),prof(:,3),'r.-',...
         prof(:,4)-50,prof(:,3),'r',...
         prof(:,4)+50,prof(:,3),'r'); hold off; grid

%%plot the (q,p) profs
for iI = 1 : 6
  iInd = (1:3) + (iI-1)*3;
  caX = caIn(iInd);
  prof = load(caX);
  semilogy(prof(:,5),prof(:,3)); hold on
  end
semilogy(prof(:,5),prof(:,3),'r.-'); hold off; grid

%%%%%%now take stuff from 70 - 120 km, in 100 layers
iUp = find(prof(:,2) >= 70);
semilogy(prof(iUp,2),prof(iUp,3),'.-'); grid
jj = iUp - iUp(1); jj = jj + 1;
%% want P(001) = prof(iUp(1),2);
%% want P(101) = prof(iUp(length(iUp)),2);
P1 = prof(iUp(1),3);
P2 = prof(iUp(length(iUp)),3);

%%%%%%%%%%%%%%%% bad %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% try j = 1 : 101, P(j=1) = 0.05, P(40) = 1.8e-3, P(101) = 2.5e-5 
%%% try Pj = Aj^2 + Bj + C
C = [1 1 1 ; 1600 40 1; 101*101 101 1];
G = [5e-2 1.8e-3 2.5e-5]';
aa = inv(C)*G;
a1 = aa(1); b1 = aa(2); c1 = aa(3);
j = 1 : 101;
PlevUp = ( a1*j.^2 + b1*j + c1 );
plot(j,PlevUp)
%%%%%%%%%%%%%%%% bad %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% good %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% try j = 1 : 101, P(j=1) = 0.05 mb, P(101) = 2.5e-5 mb
%%% try Pj = Aexp(-j/jo)
C = [1 -1 ; 1 -101];
G = [5e-2 2.5e-5]'; G = log(G);
aa = inv(C)*G;        %% -2.919723248958570e+00, 7.600902459542083e-02

a1 = exp(aa(1)); b1 = 1/aa(2);  %% 5.394861557009636e-02  1.315633249239519e+01
j = 1 : 101;
PlevUp = a1 * exp(-j/b1);
plot(j,PlevUp)

%now see how the layers etc are
hnew = interp1(log(prof(:,3)),prof(:,2),log(PlevUp),[],'extrap');
tnew = interp1(log(prof(:,3)),prof(:,4),log(PlevUp),[],'extrap');
qnew = interp1(log(prof(:,3)),prof(:,5),log(PlevUp),[],'extrap');

semilogy(prof(:,2),prof(:,3),hnew,PlevUp,'r.-'); grid on; title('Alt');
semilogy(prof(:,4),prof(:,3),tnew,PlevUp,'r.-'); grid on; title('T');
semilogy(prof(:,5),prof(:,3),qnew,PlevUp,'r.-'); grid on; title('q');

plot(hnew(1:length(hnew)-1),diff(hnew)) %%shows the layers are 0.5-0.8km thick

%%%%%%%%%%%%%%%% good %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% make the files for run7co2
px1 = PlevUp(1:100);
px2 = PlevUp(2:101);
raPavg = (px1-px2)./log(px1./px2);

raHavg = interp1(log(prof(:,3)),prof(:,2),log(raPavg),[],'extrap');
raTavg = interp1(log(prof(:,3)),prof(:,4),log(raPavg),[],'extrap');
raQavg = interp1(log(prof(:,3)),prof(:,5),log(raPavg),[],'extrap');

semilogy(prof(:,2),prof(:,3),raHavg,raPavg,'r.-'); grid on; title('Alt');
semilogy(prof(:,4),prof(:,3),raTavg,raPavg,'r.-'); grid on; title('T');
semilogy(prof(:,5),prof(:,3),raQavg,raPavg,'r.-'); grid on; title('q');

%%see SUBROUTINE read_upperatm_lte_temperature in kreadVTprofiles.f (KCARTA)
MGC = 8.314674269981136;
kAtm2mb = 1013.25;


% ppmv = (number of gas molecules)/(number of air molecules) * 1e6 
% pV = NRT ==> N(g) = p(g) V(a)/RT(g),  N(a) = p(a) V(a)/ RT(a) 
% now V(g) == V(a), T(g) == T(a) and p(total) = p(a) + p(g) 
% thus ppmv = N(g)/N(a) = p(g)/p(a) = p(g)/(p(total) - p(g)) * 1e6 
% thus p(g) = (p(total)/(1e6+ppmv)) ppmv 

press     = raPavg/kAtm2mb;                  %% in mb -> atm
partpress = press.*raQavg./(1.0e6+raQavg);   %% in atm

raLayThick = diff(hnew)*1000*100;  %% in km -> m -> cm
                                   %% NOTE use of hnew, not raHavg
co2 = raLayThick * kAtm2mb * 100.0 .* partpress;
co2 = co2./1e9/MGC./raTavg;

%%%now need to output [iLay  press  partpress raTavg co2] 
array = [1:100; press;  partpress; raTavg; co2];
fid = fopen('refgas2','w');
fprintf(fid,'%3i  %10.8e %10.8e %7.3f %10.8e \n',array);  
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% need new layering for NLTE 
semilogy(j,Plev,'b.-',j,PlevUp,'r.-')

%% combine them to make a new klayers_NLTE between 0 - 80 km
usual   = find(Plev   >= max(PlevUp));
newones = find(PlevUp >= min(Plev)); newones = newones(1:2:length(newones));

klayers_Plev = [Plev(usual) PlevUp(newones)];
heights_Plev = interp1(log(prof(:,3)),prof(:,2),log(klayers_Plev));
semilogy(heights_Plev,klayers_Plev,'.-')

%%% OR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% try 1 : 38  as the current AIRS levels (1100-300 mb)
%%    38 : 121  from 300 - 0.005 mb (84 levels)
CC = [1 -1 ; 1 -84]; 
GG = [300 5e-3 ]'; GG = log(GG); 
aaaa = inv(CC)*GG;        %% 5.808595905357528e+00, 1.048134307013276e-01
 
aa1 = exp(aaaa(1)); bb1 = 1/aaaa(2);  
  %% 3.331510215557295e+02  9.540762031247333e+00
j2 = 2 : 84; 
PlevUp2 = aa1 * exp(-j2/bb1); 
plot(j2,PlevUp2) 

klayers_Plev = [Plev(1:38) PlevUp2]; 
heights_Plev = interp1(log(prof(:,3)),prof(:,2),log(klayers_Plev)); 
semilogy(heights_Plev,klayers_Plev,'.-') 
% plot(diff(heights_Plev))
% plot(heights_Plev(1:length(heights_Plev)-1),diff(heights_Plev))

%%% OR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% try 1 : 66  as the current AIRS levels (1100-51.5 mb)
%%    66 : 121  from 51.5 - 0.005 mb (56 levels)
CC = [1 -1 ; 1 -56]; 
GG = [51.5 5e-3 ]'; GG = log(GG); 
aaaa = inv(CC)*GG;  
aa1 = exp(aaaa(1)); bb1 = 1/aaaa(2);  
j2 = 2 : 55; 
PlevUp2 = aa1 * exp(-j2/bb1); 
plot(j2,PlevUp2) 

klayers_Plev = [Plev(1:66) PlevUp2]; 
heights_Plev = interp1(log(prof(:,3)),prof(:,2),log(klayers_Plev)); 
semilogy(heights_Plev,klayers_Plev,'.-') 
% plot(diff(heights_Plev))
% plot(heights_Plev(1:length(heights_Plev)-1),diff(heights_Plev),...
%      husual(1:length(husual)-1),diff(husual))

%%% OR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% find where the current AIRS layers become 1 km thick
%%    81 : 120  from 9.5 - 0.005 mb (40 levels)
CC = [1 -1 ; 1 -40]; 
GG = [9.5 5e-3 ]'; GG = log(GG); 
aaaa = inv(CC)*GG;  
aa1 = exp(aaaa(1)); bb1 = 1/aaaa(2);  
j2 = 2 : 40; 
PlevUp2 = aa1 * exp(-j2/bb1); 
plot(j2,PlevUp2) 

klayers_Plev = [Plev(1:81) PlevUp2]; 
heights_Plev = interp1(log(prof(:,3)),prof(:,2),log(klayers_Plev)); 
semilogy(heights_Plev,klayers_Plev,'.-') 
% plot(diff(heights_Plev))
% plot(heights_Plev(1:length(heights_Plev)-1),diff(heights_Plev),...
%      husual(1:length(husual)-1),diff(husual))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% also need to make parameter files for kCARTA to know layer thickness etc
make_param_files;

