%>>>>>>>>>J-stats ATom paper#3 - Observations INPUT

% CAFS (121,123 total points)
XXlon = ncread('Jstat_obs_CAFS.nc','XXlon');
XXlat = ncread('Jstat_obs_CAFS.nc','XXlat');
XXp   = ncread('Jstat_obs_CAFS.nc','XXp');
XXsza = ncread('Jstat_obs_CAFS.nc','XXsza');
M = length(XXlon)
J1cld = ncread('Jstat_obs_CAFS.nc','J1cld');
J1clr = ncread('Jstat_obs_CAFS.nc','J1clr');
J2cld = ncread('Jstat_obs_CAFS.nc','J2cld');
J2clr = ncread('Jstat_obs_CAFS.nc','J2clr');


Jxx(1:M,1,1) = J1cld(1:M);
Jxx(1:M,1,2) = J1clr(1:M);
Jxx(1:M,2,1) = J2cld(1:M);
Jxx(1:M,2,2) = J2clr(1:M);


% now map out in the key BLOCKS:
% so far:  Block1 = Trop Pacific 20S-20N x 160E-240E
%          Block2 = NorthPacific  20N-50N x 170E-225E
% NOTE current logic will not place a block across Longitude = 0

Blocks = 3;
LAT1 = [ -20;
         +20;
         -55];
LAT2 = [ +20;
         +50
         +55];
LON1 = [ 160;
         170;
           0];
LON2 = [ 240;
         225;
         360];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SZA0 = 0.77;  %% 40 deg SZA

XJ(:,1) = J1cld./J1clr;
XJ(:,2) = J2cld./J2clr;

for j=1:2
 lnXJ = log(XJ(:,j));
 for b=1:Blocks
      lnXx = 0;
      lnXn = 500;
      Jcount0 = 0;
      Jcount(1:201) = 0;
  for i=1:M
   if XXsza(i) > SZA0
    if XXlat(i)>=LAT1(b) & XXlat(i)<=LAT2(b)
    if XXlon(i)>=LON1(b) & XXlon(i)<=LON2(b)
        Jk = min(201, max(1, 101+floor(100*lnXJ(i)+0.5)));
      Jcount0 = Jcount0+1;
      Jcount(Jk) = Jcount(Jk)+1;
      lnXx = max(lnXx,lnXJ(i));
      lnXn = min(lnXn,lnXJ(i));
    end
    end
   end
  end
    for n=1:201
     prJ(n,j,b) = Jcount(n)/Jcount0;
    end
    [j,b, Jcount0]
    [lnXx, lnXn]
 end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% blocks: b1 = 11504   B2 = 4867   B3 = 41968
%  cos(sza) > 0.77 (0-40 deg)

%%%%%%%% calculate CAFS offset from peak & plot of CAFS prob dist of ln(lcd/clr)
Xp = -1.0:0.01:1.0;

figure(101)
  hold on
  set(gcf, 'units','points','outerposition',[50 50 600 500]);
  figure(gcf)
   hold on
  title(' across J1-J2, down blocks 1-2-3-')
for j=1:2
 for b=1:Blocks
  k = j + 2*(b-1)  % index for multiple graphs
  prK = prJ(:,j,b);
  prK(1) = 0;
  dprK = prK;
  for kk=1:6
   for i=2:199
    dprK(i) = (prK(i-1)+2*prK(i)+prK(i+1))/4;
   end
    prK = dprK;
  end
   [JA,JI] = sort(prK,'descend');
   Jpeak(j,b) = JI(1);
 subplot(Blocks,2,k)
   hold on
   grid on
   set(gca,'XLim',[-.5,.5],'XTick',-.5:0.1:.5)
   plot(Xp(51:152),prJ(51:152,j,b),'.r','MarkerSize',3)
   plot(Xp(51:152),prK(51:152),'k')
   xlabel('ln(J-cldy/J-clear)')
   title('CAFS:  across J1-J2, down blocks 1-2-3-')
 end
end
Jpeak
%

offset(1) = (101-Jpeak(1,3))*0.010
offset(2) = (101-Jpeak(2,3))*0.010
%%   offset = [+0.04; 0.09]
%%
%%%%%% new with TUV better albedo and no AOD
%%  offset = [+0.0200,  +0.0100]
%%% drop the offset correction

offset = [0,0];


%%%%%%%%%%% mean CAFS Profles (CJp) of J-cldy and J-clr
%%%%% 3-blocks, high sun to avoid corrections
%%%%% Reduce J-clr (TUV) value by exp(offset)
SZA0 = 0.80

%  CJp(150:950 hPa, J1:J2, cldy:clr, block1:block3)
%       k            j      1:2         b
CJ(1:9,1:2,1:2,1:Blocks) = 0.0;
CJ2(1:9,1:2,1:2,1:Blocks)= 0.0;
countkjb(1:9,1:2,1:Blocks) = 0;

for i=1:M
 if XXsza(i) > SZA0
    k = max(1,min(9,floor(XXp(i)/100)));
  for b=1:Blocks
   if XXlat(i)>= LAT1(b) & XXlat(i)<= LAT2(b)
   if XXlon(i)>= LON1(b) & XXlon(i)<= LON2(b)
    for j=1:2
     if Jxx(i,j,2) >0 & Jxx(i,j,1) >0
      CJ(k,j,1,b) = CJ(k,j,1,b) + Jxx(i,j,1);
      CJ(k,j,2,b) = CJ(k,j,2,b) + Jxx(i,j,2);
      CJ2(k,j,1,b) = CJ2(k,j,1,b) + Jxx(i,j,1)^2;
      CJ2(k,j,2,b) = CJ2(k,j,2,b) + Jxx(i,j,2)^2;
      countkjb(k,j,b) = countkjb(k,j,b) + 1;
     end
    end
   end
   end
  end
 end
end

CJp(1:9,1:2,1:2,1:Blocks) = 0.0;
CJp2(1:9,1:2,1:2,1:Blocks)= 0.0;

for b=1:Blocks
 for cc=1:2
  for j=1:2
   for k=1:9
    if countkjb(k,j,b) > 0
  CJp(k,j,cc,b) = CJ(k,j,cc,b)./countkjb(k,j,b);
  CJp2(k,j,cc,b) = sqrt(CJ2(k,j,cc,b)./countkjb(k,j,b)-CJp(k,j,cc,b).^2);
    end
   end
  end
 end
end

% reduce the clear sky J by exp(offset(j)), not s.d.'s
for j=1:2
 for k=1:9
  for b=1:Blocks
   CJp(k,j,3,b) = CJp(k,j,2,b)*exp(-offset(j));
  end
 end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plots of avg CAFS profile J's w/ std dev

PJ = 150:100:950;
PJ(1) = NaN;

Xlb  = ['J-O1D cldy';
        'J-NO2 cldy'];
Ylb  = ['pressure (hPa)']
Xlm = [ 8e-5, 15e-3];
Llb = ['Global ';
       'N. Pac.';
       'Tr.Pac.';
       '-SD cld'
       '-SD clr'];

figure(102)
   hold on
   set(gcf, 'units','points','outerposition',[50 50 700 400]);
   figure(gcf)
   hold on
 for j=1:2
  subplot(1,2,j)
  hold on
  set(gca, 'YLim',[200, 1000], 'YTick', 200:100:1000, 'YDir', 'reverse')
  set(gca, 'XLim',[0, Xlm(j)])
  set(gca, 'FontSize',12, 'LineWidth',1)
  xlabel(Xlb(j,:))
  ylabel(Ylb)
  grid on
  box on
  plot( CJp(:,j,1,3),PJ,'k','LineWidth',3)
  plot( CJp(:,j,1,2),PJ,'r','LineWidth',3)
  plot( CJp(:,j,1,1),PJ,'b','LineWidth',3)
  plot( CJp2(:,j,1,1),PJ,'*b','LineWidth',2,'MarkerSize',09)
  plot( CJp2(:,j,2,1),PJ,'ob','LineWidth',2,'MarkerSize',09)
  legend(Llb(1:5,:))
  legend('Location','North')
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%% Generate ln(cldy/clr) PDFs for 100 hPa levels
%%%%%% offset included %%%%
%  pdJ (1:pdbins, 150:950 hPa, J1:J2, block1:block2)
%  pdJ ( 1:201,    1:9,         1:2,    1:2ormore)
%  Jxx (1:#CAFS, J1:J2, cldy:clr)
%  Jxx (1:M,      1:2,    1:2)

pdJ(1:201,1:9,1:2,1:Blocks) = 0;

SZA0 = 0.5;
%%%% this wider cos(SZA) limit covers most of ATom-1 CAFS
cpwt = 1;

for i=1:M
 if XXsza(i) > SZA0
    k = min(9,floor(XXp(i)/100));
  for b=1:Blocks
   if XXlat(i)>= LAT1(b) & XXlat(i)<= LAT2(b)
   if XXlon(i)>= LON1(b) & XXlon(i)<= LON2(b)
    for j=1:2
     if Jxx(i,j,2) >0  & Jxx(i,j,1) >0
         lnJ = log(Jxx(i,j,1)/Jxx(i,j,2)) + offset(j);
         nJ = min(201, max(1, 101+floor(100*lnJ+0.5)));
        pdJ(nJ,k,j,b) = pdJ(nJ,k,j,b) + cpwt;
     end
    end
   end
   end
  end
 end
end

%%%%%%%%%%%%%%%%%re-normalize, keep weighting for each p-level
for k=1:9
 for b=1:Blocks
  for j=1:2
    pdW(k,j,b) = sum(pdJ(:,k,j,b));
   if pdW(k,j,b) < 1
    pdW(k,j,b) = 1;
   end
   for n=1:201
    pdJ(n,k,j,b) = pdJ(n,k,j,b)/pdW(k,j,b);
   end
%%% calculate fraction < -0.50 and > 0.50 (not seen in plot)
     pdJtop(k,j,b) = sum(pdJ(152:201,k,j,b));
     pdJbot(k,j,b) = sum(pdJ(  1:50 ,k,j,b));
  end
 end
end

%%%%%%%%%%%%%%%%%%% pdJs(mooth) = smoothed binning  6 x 1-2-1
     p121(1:201)=0;
for b=1:Blocks
 for j=1:2
  for k=1:9
   for n=1:201
     p121(n) = pdJ(n,k,j,b);
   end
     p121(1) = 0;
   for m=1:6
    for n=2:200
     p121(n) = (p121(n-1)+2*p121(n)+p121(n+1))/4;
    end
   end
     pdJs(:,k,j,b) = p121(:);
  end
 end
end

%%%%% generate 3 pressure regions and averages
Plev1 = [1;
         4;
         9];
Plev2 = [3;
         8;
         9];
Plbl =['100-300';
       '300-900';
       '900-srf'];

JYstat(1:201,1:3,1:2,1:3)=0;  % (n, k3, j, b)
for k3=1:3
 for j=1:2
  for b=1:Blocks
    sumJ(1:201) = 0;
    cntJ = 0;
   for k = Plev1(k3):Plev2(k3)
    for n=1:201
     sumJ(n) = sumJ(n) + pdJs(n,k,j,b)*pdW(k,j,b);
    end
     cntJ = cntJ + pdW(k,j,b);
   end
    for n=1:201
     JYstat(n,k3,j,b) = sumJ(n)/cntJ;
    end
  end
 end
end

Jtitle = ['J-O1D Tr.Pac';
          'J-NO2 Tr.Pac';
          'J-O1D N. Pac';
          'J-NO2 N. Pac';
          'J-O1D Global';
          'J-NO2 Global'];



%%%%%%%%%%%% ln cld2clr, offset & smoothed: pdJ => pdJs => JYstat
figure(103)
  hold on
  set(gcf, 'units','points','outerposition',[50 50 800 600]);
  figure(gcf)
   hold on
for j=1:2
 for b=1:Blocks
  k = j + 2*(b-1); % index for multiple graphs
 subplot(Blocks,2,k)
   set(gca, 'FontSize',12, 'LineWidth',1)
   hold on
   box on
   grid on
   set(gca,'XLim',[-.5,.5],'XTick',-.5:0.1:.5)
   set(gca,'YLim',[0, 0.10])
   plot(Xp(51:152),JYstat(51:152,1,j,b),'b','LineWidth',2)
   plot(Xp(51:152),JYstat(51:152,2,j,b),'r','LineWidth',2)
   plot(Xp(51:152),JYstat(51:152,3,j,b),'k','LineWidth',2)
   if b==3
    xlabel('ln(J-cloudy/J-clear)')
   end
   title(Jtitle(k,:))
   legend({Plbl(1:3,:)},'FontSize',11,'Location','NorthWest')
 end
end
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%% ln cld2clr, not offset, but smoothed: prJ & prK
figure(104)
  hold on
  set(gcf, 'units','points','outerposition',[50 50 800 600]);
  figure(gcf)
   hold on
for j=1:2
 for b=1:Blocks
  k = j + 2*(b-1)  % index for multiple graphs
  prK = prJ(:,j,b);
  prK(1) = 0;
  dprK = prK;
  for kk=1:6
   for i=2:199
    dprK(i) = (prK(i-1)+2*prK(i)+prK(i+1))/4;
   end
    prK = dprK;
  end
 subplot(Blocks,2,k)
   set(gca, 'FontSize',12, 'LineWidth',1)
   hold on
   box on
   grid on
   set(gca,'XLim',[-.5,.5],'XTick',-.5:0.1:.5)
   set(gca,'YLim',[0, 0.10])
   plot(Xp(51:152),prJ(51:152,j,b),'.r','MarkerSize',5,'LineWidth',2)
   plot(Xp(51:152),prK(51:152),'k','LineWidth',2)
   if b==3
    xlabel('ln(J-cloudy/J-clear)')
   end
   title(Jtitle(k,:))
 end
end
hold off
%