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

Mdls = 10;

%model	        I     J    K
%  A   GC     144    91   37
%  B   GISS   144    90   28
%  C   UCI    320   160   34
%  D   NCAR   576   384   36
%  E   UKCA   192   144   42
%  F   MOCA   360   180   33
%  G   IFS    512   256   82
%  H   GFDL   720   360   25
%  I   UCIb   320   160   34    Briegleb-averaged clouds
%  K   GMI    288   181   37

% note that model levels P < 70 hPa are truncated
% models with reverse layers (#1 = top) have been inverted

mdlname = ['GC  ';
           'GFDL';
           'GISS';
           'GMI ';
           'IFS ';
           'MOCA';
           'NCAR';
           'UCI ';
           'UKCA';
           'UCIb'];

% GC (GEOS-Chem)
XlonA = ncread('Jstat_Amodel_GC.nc','lon');
YlatA = ncread('Jstat_Amodel_GC.nc','lat');
PlevAx= ncread('Jstat_Amodel_GC.nc','lev');
NXA = length(XlonA)
NYA = length(YlatA)
NPA = length(PlevAx)
PlevA = ncread('Jstat_Amodel_GC.nc','PlevA');
J1A = ncread('Jstat_Amodel_GC.nc','J1A');
J1Ac= ncread('Jstat_Amodel_GC.nc','J1Ac');
J2A = ncread('Jstat_Amodel_GC.nc','J2A');
J2Ac= ncread('Jstat_Amodel_GC.nc','J2Ac');

% GISS
XlonB = ncread('Jstat_Bmodel_GISS.nc','lon');
YlatB = ncread('Jstat_Bmodel_GISS.nc','lat');
PlevBx= ncread('Jstat_Bmodel_GISS.nc','lev');
NXB = length(XlonB)
NYB = length(YlatB)
NPB = length(PlevBx)
PlevB = ncread('Jstat_Bmodel_GISS.nc','PlevB');
J1B = ncread('Jstat_Bmodel_GISS.nc','J1B');
J1Bc= ncread('Jstat_Bmodel_GISS.nc','J1Bc');
J2B = ncread('Jstat_Bmodel_GISS.nc','J2B');
J2Bc= ncread('Jstat_Bmodel_GISS.nc','J2Bc');

% UCI
XlonC = ncread('Jstat_Cmodel_UCI.nc','lon');
YlatC = ncread('Jstat_Cmodel_UCI.nc','lat');
PlevCx= ncread('Jstat_Cmodel_UCI.nc','lev');
NXC = length(XlonC)
NYC = length(YlatC)
NPC = length(PlevCx)
PlevC = ncread('Jstat_Cmodel_UCI.nc','PlevC');
J1C = ncread('Jstat_Cmodel_UCI.nc','J1C');
J1Cc= ncread('Jstat_Cmodel_UCI.nc','J1Cc');
J2C = ncread('Jstat_Cmodel_UCI.nc','J2C');
J2Cc= ncread('Jstat_Cmodel_UCI.nc','J2Cc');

% NCAR
XlonD = ncread('Jstat_Dmodel_NCAR.nc','lon');
YlatD = ncread('Jstat_Dmodel_NCAR.nc','lat');
PlevDx= ncread('Jstat_Dmodel_NCAR.nc','lev');
NXD = length(XlonD)
NYD = length(YlatD)
NPD = length(PlevDx)
PlevD = ncread('Jstat_Dmodel_NCAR.nc','PlevD');
J1D = ncread('Jstat_Dmodel_NCAR.nc','J1D');
J1Dc= ncread('Jstat_Dmodel_NCAR.nc','J1Dc');
J2D = ncread('Jstat_Dmodel_NCAR.nc','J2D');
J2Dc= ncread('Jstat_Dmodel_NCAR.nc','J2Dc');

% UKCA
XlonE = ncread('Jstat_Emodel_UKCA.nc','lon');
YlatE = ncread('Jstat_Emodel_UKCA.nc','lat');
PlevEx= ncread('Jstat_Emodel_UKCA.nc','lev');
NXE = length(XlonE)
NYE = length(YlatE)
NPE = length(PlevEx)
PlevE = ncread('Jstat_Emodel_UKCA.nc','PlevE');
J1E = ncread('Jstat_Emodel_UKCA.nc','J1E');
J1Ec= ncread('Jstat_Emodel_UKCA.nc','J1Ec');
J2E = ncread('Jstat_Emodel_UKCA.nc','J2E');
J2Ec= ncread('Jstat_Emodel_UKCA.nc','J2Ec');

% MOCA
XlonF = ncread('Jstat_Fmodel_MOCA.nc','lon');
YlatF = ncread('Jstat_Fmodel_MOCA.nc','lat');
PlevFx= ncread('Jstat_Fmodel_MOCA.nc','lev');
NXF = length(XlonF)
NYF = length(YlatF)
NPF = length(PlevFx)
PlevF = ncread('Jstat_Fmodel_MOCA.nc','PlevF');
J1F = ncread('Jstat_Fmodel_MOCA.nc','J1F');
J1Fc= ncread('Jstat_Fmodel_MOCA.nc','J1Fc');
J2F = ncread('Jstat_Fmodel_MOCA.nc','J2F');
J2Fc= ncread('Jstat_Fmodel_MOCA.nc','J2Fc');

% IFS truncate levels; reverse from 1:137 to 1:82; reverse latitude
XlonG = ncread('Jstat_Gmodel_IFS.nc','lon');
YlatG = ncread('Jstat_Gmodel_IFS.nc','lat');
PlevGx= ncread('Jstat_Gmodel_IFS.nc','lev');
NXG = length(XlonG)
NYG = length(YlatG)
NPG = length(PlevGx)
PlevG = ncread('Jstat_Gmodel_IFS.nc','PlevG');
J1G = ncread('Jstat_Gmodel_IFS.nc','J1G');
J1Gc= ncread('Jstat_Gmodel_IFS.nc','J1Gc');
J2G = ncread('Jstat_Gmodel_IFS.nc','J2G');
J2Gc= ncread('Jstat_Gmodel_IFS.nc','J2Gc');

% GFDL
XlonH = ncread('Jstat_Hmodel_GFDL.nc','lon');
YlatH = ncread('Jstat_Hmodel_GFDL.nc','lat');
PlevHx= ncread('Jstat_Hmodel_GFDL.nc','lev');
NXH = length(XlonH)
NYH = length(YlatH)
NPH = length(PlevHx)
PlevH = ncread('Jstat_Hmodel_GFDL.nc','PlevH');
J1H = ncread('Jstat_Hmodel_GFDL.nc','J1H');
J1Hc= ncread('Jstat_Hmodel_GFDL.nc','J1Hc');
J2H = ncread('Jstat_Hmodel_GFDL.nc','J2H');
J2Hc= ncread('Jstat_Hmodel_GFDL.nc','J2Hc');

% UCIb = Briegleb-averaged clouds
XlonI = ncread('Jstat_Imodel_UCIb.nc','lon');
YlatI = ncread('Jstat_Imodel_UCIb.nc','lat');
PlevIx= ncread('Jstat_Imodel_UCIb.nc','lev');
NXI = length(XlonI)
NYI = length(YlatI)
NPI = length(PlevIx)
PlevI = ncread('Jstat_Imodel_UCIb.nc','PlevI');
J1I = ncread('Jstat_Imodel_UCIb.nc','J1I');
J1Ic= ncread('Jstat_Imodel_UCIb.nc','J1Ic');
J2I = ncread('Jstat_Imodel_UCIb.nc','J2I');
J2Ic= ncread('Jstat_Imodel_UCIb.nc','J2Ic');

% GMI
XlonK = ncread('Jstat_Kmodel_GMI.nc','lon');
YlatK = ncread('Jstat_Kmodel_GMI.nc','lat');
PlevKx= ncread('Jstat_Kmodel_GMI.nc','lev');
NXK = length(XlonK)
NYK = length(YlatK)
NPK = length(PlevKx)
PlevK = ncread('Jstat_Kmodel_GMI.nc','PlevK');
J1K = ncread('Jstat_Kmodel_GMI.nc','J1K');
J1Kc= ncread('Jstat_Kmodel_GMI.nc','J1Kc');
J2K = ncread('Jstat_Kmodel_GMI.nc','J2K');
J2Kc= ncread('Jstat_Kmodel_GMI.nc','J2Kc');

% pre-calculated cos(SZA) for 1x1 grid (360,180) for mid-August
cosSZA = ncread('Jstat_cosSZA.nc','cosSZA');


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

%>>>>>>>>>Process arrays etc. define weights and boundaries
% collect mean P profiles only for: (60S-60N, all Long, all t's)
    N6 = floor(NYA/6 + 1);
  PlevAx = (squeeze(mean(mean(PlevA(:,N6:(NYA-N6),:),1),2)));
    N6 = floor(NYB/6 + 1);
  PlevBx = (squeeze(mean(mean(PlevB(:,N6:(NYB-N6),:),1),2)));
    N6 = floor(NYC/6 + 1);
  PlevCx = (squeeze(mean(mean(PlevC(:,N6:(NYC-N6),:),1),2)));
    N6 = floor(NYD/6 + 1);
  PlevDx = (squeeze(mean(mean(PlevD(:,N6:(NYD-N6),:),1),2)));
    N6 = floor(NYE/6 + 1);
  PlevEx = (squeeze(mean(mean(PlevE(:,N6:(NYE-N6),:),1),2)));
    N6 = floor(NYF/6 + 1);
  PlevFx = (squeeze(mean(mean(PlevF(:,N6:(NYF-N6),:),1),2)));
    N6 = floor(NYG/6 + 1);
  PlevGx = (squeeze(mean(mean(PlevG(:,N6:(NYG-N6),:),1),2)));
    N6 = floor(NYH/6 + 1);
  PlevHx = (squeeze(mean(mean(PlevH(:,N6:(NYH-N6),:),1),2)));
    N6 = floor(NYI/6 + 1);
  PlevIx = (squeeze(mean(mean(PlevI(:,N6:(NYI-N6),:),1),2)));
    N6 = floor(NYK/6 + 1);
  PlevKx = (squeeze(mean(mean(PlevK(:,N6:(NYK-N6),:),1),2)));

% set up weightings for each level (delta-P in hPa, scaled out later)
 for l=2:NPA-1
  pwtA(l) = 0.5*(PlevAx(l-1) - PlevAx(l+1));
 end
  pwtA(1) = pwtA(2)      ;
  pwtA(NPA) = pwtA(NPA-1);
 for l=2:NPB-1
  pwtB(l) = 0.5*(PlevBx(l-1) - PlevBx(l+1));
 end
  pwtB(1) = pwtB(2)      ;
  pwtB(NPB) = pwtB(NPB-1);
 for l=2:NPC-1
  pwtC(l) = 0.5*(PlevCx(l-1) - PlevCx(l+1));
 end
  pwtC(1) = pwtC(2)      ;
  pwtC(NPC) = pwtC(NPC-1);
 for l=2:NPD-1
  pwtD(l) = 0.5*(PlevDx(l-1) - PlevDx(l+1));
 end
  pwtD(1) = pwtD(2)      ;
  pwtD(NPD) = pwtD(NPD-1);
 for l=2:NPE-1
  pwtE(l) = 0.5*(PlevEx(l-1) - PlevEx(l+1));
 end
  pwtE(1) = pwtE(2)      ;
  pwtE(NPE) = pwtE(NPE-1);
 for l=2:NPF-1
  pwtF(l) = 0.5*(PlevFx(l-1) - PlevFx(l+1));
 end
  pwtF(1) = pwtF(2)      ;
  pwtF(NPF) = pwtF(NPF-1);
 for l=2:NPG-1
  pwtG(l) = 0.5*(PlevGx(l-1) - PlevGx(l+1));
 end
  pwtG(1) = pwtG(2)      ;
  pwtG(NPG) = pwtG(NPG-1);
 for l=2:NPH-1
  pwtH(l) = 0.5*(PlevHx(l-1) - PlevHx(l+1));
 end
  pwtH(1) = pwtH(2)      ;
  pwtH(NPH) = pwtH(NPH-1);
 for l=2:NPI-1
  pwtI(l) = 0.5*(PlevIx(l-1) - PlevIx(l+1));
 end
  pwtI(1) = pwtI(2)      ;
  pwtI(NPI) = pwtI(NPI-1);
 for l=2:NPK-1
  pwtK(l) = 0.5*(PlevKx(l-1) - PlevKx(l+1));
 end
  pwtK(1) = pwtK(2)      ;
  pwtK(NPK) = pwtK(NPK-1);

% define cosine(latitude) weights
  cwtA = cos(YlatA(:)*pi/180.);
  cwtB = cos(YlatB(:)*pi/180.);
  cwtC = cos(YlatC(:)*pi/180.);
  cwtD = cos(YlatD(:)*pi/180.);
  cwtE = cos(YlatE(:)*pi/180.);
  cwtF = cos(YlatF(:)*pi/180.);
  cwtG = cos(YlatG(:)*pi/180.);
  cwtH = cos(YlatH(:)*pi/180.);
  cwtI = cos(YlatI(:)*pi/180.);
  cwtK = cos(YlatK(:)*pi/180.);

% locate solar zenith angles for each i,j and t(=1:24h) UT
%%% assumes all data are in form 24 hrs of J's in Aug
%%% pre-calc cos(SZA) for 1x1 grid (360,180)

% make a lat-lng grid for each model grid giving  1x1-deg cell number
% then make a cosSZA grid at UTT= 00h for this grid
for j=1:NYA
  jj = floor(YlatA(j) + 90.5);
  j1x1A(j) = min(180,max(1,jj));
end
for i=1:NXA
  ii = floor(XlonA(i) + 0.5) ;
  i1x1A(i) = min(360,max(1,ii));
end

for j=1:NYB
  jj = floor(YlatB(j) + 90.5);
  j1x1B(j) = min(180,max(1,jj));
end
for i=1:NXB
  ii = floor(XlonB(i) + 0.5) ;
  i1x1B(i) = min(360,max(1,ii));
end

for j=1:NYC
  jj = floor(YlatC(j) + 90.5);
  j1x1C(j) = min(180,max(1,jj));
end
for i=1:NXC
  ii = floor(XlonC(i) + 0.5) ;
  i1x1C(i) = min(360,max(1,ii));
end

for j=1:NYD
  jj = floor(YlatD(j) + 90.5);
  j1x1D(j) = min(180,max(1,jj));
end
for i=1:NXD
  ii = floor(XlonD(i) + 0.5) ;
  i1x1D(i) = min(360,max(1,ii));
end

for j=1:NYE
  jj = floor(YlatE(j) + 90.5);
  j1x1E(j) = min(180,max(1,jj));
end
for i=1:NXE
  ii = floor(XlonE(i) + 0.5) ;
  i1x1E(i) = min(360,max(1,ii));
end

for j=1:NYF
  jj = floor(YlatF(j) + 90.5);
  j1x1F(j) = min(180,max(1,jj));
end
for i=1:NXF
  ii = floor(XlonF(i) + 0.5) ;
  i1x1F(i) = min(360,max(1,ii));
end

for j=1:NYG
  jj = floor(YlatG(j) + 90.5);
  j1x1G(j) = min(180,max(1,jj));
end
for i=1:NXG
  ii = floor(XlonG(i) + 0.5) ;
  i1x1G(i) = min(360,max(1,ii));
end

for j=1:NYH
  jj = floor(YlatH(j) + 90.5);
  j1x1H(j) = min(180,max(1,jj));
end
for i=1:NXH
  ii = floor(XlonH(i) + 0.5) ;
  i1x1H(i) = min(360,max(1,ii));
end

for j=1:NYI
  jj = floor(YlatI(j) + 90.5);
  j1x1I(j) = min(180,max(1,jj));
end
for i=1:NXI
  ii = floor(XlonI(i) + 0.5) ;
  i1x1I(i) = min(360,max(1,ii));
end

for j=1:NYK
  jj = floor(YlatK(j) + 90.5);
  j1x1K(j) = min(180,max(1,jj));
end
for i=1:NXK
  ii = floor(XlonK(i) + 0.5) ;
  i1x1K(i) = min(360,max(1,ii));
end


for j=1:NYA
 for i=1:NXA
  cosA(i,j) = cosSZA(i1x1A(i),j1x1A(j));
 end
end
for j=1:NYB
 for i=1:NXB
  cosB(i,j) = cosSZA(i1x1B(i),j1x1B(j));
 end
end
for j=1:NYC
 for i=1:NXC
  cosC(i,j) = cosSZA(i1x1C(i),j1x1C(j));
 end
end
for j=1:NYD
 for i=1:NXD
  cosD(i,j) = cosSZA(i1x1D(i),j1x1D(j));
 end
end
for j=1:NYE
 for i=1:NXE
  cosE(i,j) = cosSZA(i1x1E(i),j1x1E(j));
 end
end
for j=1:NYF
 for i=1:NXF
  cosF(i,j) = cosSZA(i1x1F(i),j1x1F(j));
 end
end
for j=1:NYG
 for i=1:NXG
  cosG(i,j) = cosSZA(i1x1G(i),j1x1G(j));
 end
end
for j=1:NYH
 for i=1:NXH
  cosH(i,j) = cosSZA(i1x1H(i),j1x1H(j));
 end
end
for j=1:NYI
 for i=1:NXI
  cosI(i,j) = cosSZA(i1x1I(i),j1x1I(j));
 end
end
for j=1:NYK
 for i=1:NXK
  cosK(i,j) = cosSZA(i1x1K(i),j1x1K(j));
 end
end


% now map out (i,j) pixels 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];

blockA(1:NXA,1:NYA,1:Blocks) = int16(0);
for j=1:NYA
 for i=1:NXA
  for b=1:Blocks
   if YlatA(j)>=LAT1(b) & YlatA(j)<=LAT2(b)
    if XlonA(i)>=LON1(b) & XlonA(i)<=LON2(b)
      blockA(i,j,b) = 1;
    end
   end
  end
 end
end

blockB(1:NXB,1:NYB,1:Blocks) = int16(0);
for j=1:NYB
 for i=1:NXB
  for b=1:Blocks
   if YlatB(j)>=LAT1(b) & YlatB(j)<=LAT2(b)
    if XlonB(i)>=LON1(b) & XlonB(i)<=LON2(b)
      blockB(i,j,b) = 1;
    end
   end
  end
 end
end

blockC(1:NXC,1:NYC,1:Blocks) = int16(0);
for j=1:NYC
 for i=1:NXC
  for b=1:Blocks
   if YlatC(j)>=LAT1(b) & YlatC(j)<=LAT2(b)
    if XlonC(i)>=LON1(b) & XlonC(i)<=LON2(b)
      blockC(i,j,b) = 1;
    end
   end
  end
 end
end

blockD(1:NXD,1:NYD,1:Blocks) = int16(0);
for j=1:NYD
 for i=1:NXD
  for b=1:Blocks
   if YlatD(j)>=LAT1(b) & YlatD(j)<=LAT2(b)
    if XlonD(i)>=LON1(b) & XlonD(i)<=LON2(b)
      blockD(i,j,b) = 1;
    end
   end
  end
 end
end

blockE(1:NXE,1:NYE,1:Blocks) = int16(0);
for j=1:NYE
 for i=1:NXE
  for b=1:Blocks
   if YlatE(j)>=LAT1(b) & YlatE(j)<=LAT2(b)
    if XlonE(i)>=LON1(b) & XlonE(i)<=LON2(b)
      blockE(i,j,b) = 1;
    end
   end
  end
 end
end

blockF(1:NXF,1:NYF,1:Blocks) = int16(0);
for j=1:NYF
 for i=1:NXF
  for b=1:Blocks
   if YlatF(j)>=LAT1(b) & YlatF(j)<=LAT2(b)
    if XlonF(i)>=LON1(b) & XlonF(i)<=LON2(b)
      blockF(i,j,b) = 1;
    end
   end
  end
 end
end

blockG(1:NXG,1:NYG,1:Blocks) = int16(0);
for j=1:NYG
 for i=1:NXG
  for b=1:Blocks
   if YlatG(j)>=LAT1(b) & YlatG(j)<=LAT2(b)
    if XlonG(i)>=LON1(b) & XlonG(i)<=LON2(b)
      blockG(i,j,b) = 1;
    end
   end
  end
 end
end

blockH(1:NXH,1:NYH,1:Blocks) = int16(0);
for j=1:NYH
 for i=1:NXH
  for b=1:Blocks
   if YlatH(j)>=LAT1(b) & YlatH(j)<=LAT2(b)
    if XlonH(i)>=LON1(b) & XlonH(i)<=LON2(b)
      blockH(i,j,b) = 1;
    end
   end
  end
 end
end

blockI(1:NXI,1:NYI,1:Blocks) = int16(0);
for j=1:NYI
 for i=1:NXI
  for b=1:Blocks
   if YlatI(j)>=LAT1(b) & YlatI(j)<=LAT2(b)
    if XlonI(i)>=LON1(b) & XlonI(i)<=LON2(b)
      blockI(i,j,b) = 1;
    end
   end
  end
 end
end

blockK(1:NXK,1:NYK,1:Blocks) = int16(0);
for j=1:NYK
 for i=1:NXK
  for b=1:Blocks
   if YlatK(j)>=LAT1(b) & YlatK(j)<=LAT2(b)
    if XlonK(i)>=LON1(b) & XlonK(i)<=LON2(b)
      blockK(i,j,b) = 1;
    end
   end
  end
 end
end


%%%% use std 100 hPa intervals as pressure blocks, add together later
for l=1:NPA
   k = min(9,max(floor(PlevAx(l)/100),1));
   blockpA(l) = k;
end
for l=1:NPB
   k = min(9,max(floor(PlevBx(l)/100),1));
   blockpB(l) = k;
end
for l=1:NPC
   k = min(9,max(floor(PlevCx(l)/100),1));
   blockpC(l) = k;
end
for l=1:NPD
   k = min(9,max(floor(PlevDx(l)/100),1));
   blockpD(l) = k;
end
for l=1:NPE
   k = min(9,max(floor(PlevEx(l)/100),1));
   blockpE(l) = k;
end
for l=1:NPF
   k = min(9,max(floor(PlevFx(l)/100),1));
   blockpF(l) = k;
end
for l=1:NPG
   k = min(9,max(floor(PlevGx(l)/100),1));
   blockpG(l) = k;
end
for l=1:NPH
   k = min(9,max(floor(PlevHx(l)/100),1));
   blockpH(l) = k;
end
for l=1:NPI
   k = min(9,max(floor(PlevIx(l)/100),1));
   blockpI(l) = k;
end
for l=1:NPK
   k = min(9,max(floor(PlevKx(l)/100),1));
   blockpK(l) = k;
end


figure(200)
contour(cosSZA')
colorbar

%%%% need to reorder model # based on alphabetic name
%  01 A GC     144    91    37
%  02 H GFDL   720   360    25
%  03 B GISS   144    90    28
%  04 K GMI    288   181    37
%  05 G IFS    512   256    82
%  06 F MOCA   360   180    33
%  07 D NCAR   576   384    36
%  08 C UCI    320   160    34
%  09 E UKCA   192   144    42
%  10 I UCIb   320   160    34

mdlname = ['GC  ';
           'GFDL';
           'GISS';
           'GMI ';
           'IFS ';
           'MOCA';
           'NCAR';
           'UCI ';
           'UKCA';
           'UCIb'];


for m=1:Mdls

if m==1
  cosX = cosA;
  blockX = blockA;
  J1 = J1Ac(:,:,1,24);
  J2 = J2Ac(:,:,1,24);
end
if m==2
  cosX = cosH;
  blockX = blockH;
  J1 = J1Hc(:,:,1,24);
  J2 = J2Hc(:,:,1,24);
end
if m==3
  cosX = cosB;
  blockX = blockB;
  J1 = J1Bc(:,:,1,24);
  J2 = J2Bc(:,:,1,24);
end
if m==4
  cosX = cosK;
  blockX = blockK;
  J1 = J1Kc(:,:,1,24);
  J2 = J2Kc(:,:,1,24);
end
if m==5
  cosX = cosG;
  blockX = blockG;
  J1 = J1Gc(:,:,1,24);
  J2 = J2Gc(:,:,1,24);
end
if m==6
  cosX = cosF;
  blockX = blockF;
  J1 = J1Fc(:,:,1,24);
  J2 = J2Fc(:,:,1,24);
end
if m==7
  cosX = cosD;
  blockX = blockD;
  J1 = J1Dc(:,:,1,24);
  J2 = J2Dc(:,:,1,24);
end
if m==8
  cosX = cosC;
  blockX = blockC;
  J1 = J1Cc(:,:,1,24);
  J2 = J2Cc(:,:,1,24);
end
if m==9
  cosX = cosE;
  blockX = blockE;
  J1 = J1Ec(:,:,1,24);
  J2 = J2Ec(:,:,1,24);
end
if m==10
  cosX = cosI;
  blockX = blockI;
  J1 = J1Ic(:,:,1,24);
  J2 = J2Ic(:,:,1,24);
end


figure(200+m)
    set(gcf, 'units','points','outerposition',[200 200 1000 200]);
    subplot(1,4,1)
        title(mdlname(m,:),'FontSize',16)
        xlabel('cos(SZA)')
        hold on
        contour(cosX')
        colorbar
    subplot(1,4,2)
        contour(J2')
        colorbar
        xlabel('J-NO2 at 0 km (clear)')
    subplot(1,4,3)
        contour(J1')
        colorbar
        xlabel('J-O1D at 0 km (clear)')
    subplot(1,4,4)
        hold on
        contour(blockX(:,:,1)')
        contour(blockX(:,:,2)')
        colorbar
        xlabel('Trop Pac and North Pac blocks')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%