%>>>>>>>>J-stats ATom paper#3 - Generate model J stats
%   SLOW  SLOW  SLOW

%%%% need to reorder model # based on 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 switched to Briegleb-averaged = 'b'

Mdls = 10

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

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

% calc model avg J's for cosSZA > 0.8 to match CAFS
Jprof(1:99,1:2,1:2,1:Blocks,1:Mdls)=0;
Pprof(1:99,1:Blocks,1:Mdls)=0;

pdMJ(1:201,1:9,1:2,1:Blocks,1:Mdls) = 0;
JXstat(1:201,1:3,1:2,1:Blocks,1:Mdls)=0;


SZA0 = 0.80

%%%%%%
for model = 1:Mdls
%%%%%%

 if model == 1
   NX = NXA;
   NY = NYA;
   NP = NPA;
  Mlb(model,:),[model,NX,NY,NP]
   block = blockA;
   blockp = blockpA;
   plevel = PlevA;
   cosy = cosA;
   pwt = pwtA;
   cwt = cwtA;
   J11 = J1A;
   J12 = J1Ac;
   J21 = J2A;
   J22 = J2Ac;
 end
 if model == 3
   NX = NXB;
   NY = NYB;
   NP = NPB;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockB;
   blockp = blockpB;
   plevel = PlevB;
   cosy = cosB;
   pwt = pwtB;
   cwt = cwtB;
   J11 = J1B;
   J12 = J1Bc;
   J21 = J2B;
   J22 = J2Bc;
 end
 if model == 7
   NX = NXD;
   NY = NYD;
   NP = NPD;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockD;
   blockp = blockpD;
   plevel = PlevD;
   cosy = cosD;
   pwt = pwtD;
   cwt = cwtD;
   J11 = J1D;
   J12 = J1Dc;
   J21 = J2D;
   J22 = J2Dc;
 end
 if model == 8
  NX = NXC;
  NY = NYC;
  NP = NPC;
  Mlb(model,:),[model,NX,NY,NP]
  block = blockC;
  blockp = blockpC;
  plevel = PlevC;
  cosy = cosC;
  pwt = pwtC;
  cwt = cwtC;
  J11 = J1C;
  J12 = J1Cc;
  J21 = J2C;
  J22 = J2Cc;
 end
 if model == 9
   NX = NXE;
   NY = NYE;
   NP = NPE;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockE;
   blockp = blockpE;
   plevel = PlevE;
   cosy = cosE;
   pwt = pwtE;
   cwt = cwtE;
   J11 = J1E;
   J12 = J1Ec;
   J21 = J2E;
   J22 = J2Ec;
 end
 if model == 6
   NX = NXF;
   NY = NYF;
   NP = NPF
   Mlb(model,:),[model,NX,NY,NP]
   block = blockF;
   blockp = blockpF;
   plevel = PlevF;
   cosy = cosF;
   pwt = pwtF;
   cwt = cwtF;
   J11 = J1F;
   J12 = J1Fc;
   J21 = J2F;
   J22 = J2Fc;
 end
 if model == 5
   NX = NXG;
   NY = NYG;
   NP = NPG;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockG;
   blockp = blockpG;
   plevel = PlevG;
   cosy = cosG;
   pwt = pwtG;
   cwt = cwtG;
   J11 = J1G;
   J12 = J1Gc;
   J21 = J2G;
   J22 = J2Gc;
 end
 if model == 2
   NX = NXH;
   NY = NYH;
   NP = NPH;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockH;
   blockp = blockpH;
   plevel = PlevH;
   cosy = cosH;
   pwt = pwtH;
   cwt = cwtH;
   J11 = J1H;
   J12 = J1Hc;
   J21 = J2H;
   J22 = J2Hc;
 end
 if model == 10
   NX = NXI;
   NY = NYI;
   NP = NPI;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockI;
   blockp = blockpI;
   plevel = PlevI;
   cosy = cosI;
   pwt = pwtI;
   cwt = cwtI;
   J11 = J1I;
   J12 = J1Ic;
   J21 = J2I;
   J22 = J2Ic;
 end
 if model == 4
   NX = NXK;
   NY = NYK;
   NP = NPK;
   Mlb(model,:),[model,NX,NY,NP]
   block = blockK;
   blockp = blockpK;
   plevel = PlevK;
   cosy = cosK;
   pwt = pwtK;
   cwt = cwtK;
   J11 = J1K;
   J12 = J1Kc;
   J21 = J2K;
   J22 = J2Kc;
 end

%%%%%% core code, same for all models
 J11b(1:90,1:Blocks)=0; % J1 cldy
 J12b(1:90,1:Blocks)=0; % J1 clear
 J21b(1:90,1:Blocks)=0; % J2 cldy
 J22b(1:90,1:Blocks)=0;
 Pb(1:90,1:Blocks)=0;
 Nb(1:Blocks)=0;
 for i=1:NX
  for j=1:NY
   for b=1:Blocks
    if block(i,j,b) == 1
     for t=1:24
        ishift = floor(t*NX/24.);
        ii = i + ishift;
      if ii>NX
        ii = ii-NX;
      end
      if cosy(ii,j) > SZA0
        for l=1:NP
         J11b(l,b) = J11b(l,b) + J11(i,j,l,t);
         J12b(l,b) = J12b(l,b) + J12(i,j,l,t);
         J21b(l,b) = J21b(l,b) + J21(i,j,l,t);
         J22b(l,b) = J22b(l,b) + J22(i,j,l,t);
         Pb(l,b) = Pb(l,b) + plevel(i,j,l);
        end
         Nb(b) = Nb(b) + 1;
       end
     end
    end
   end
  end
 end

% JProf (model-level,J#,cld:clr,Block,model#) (k,j,cc,b)
 for b=1:Blocks
  for k=1:NP
   Jprof(k,1,1,b,model) = J11b(k,b)./Nb(b);
   Jprof(k,1,2,b,model) = J12b(k,b)./Nb(b);
   Jprof(k,2,1,b,model) = J21b(k,b)./Nb(b);
   Jprof(k,2,2,b,model) = J22b(k,b)./Nb(b);
   Pprof(k,b,model) = Pb(k,b)./Nb(b);
  end
 end

%%%%%%%% set up pd's of Jcldy/Jclear %%%%%%%%%%%

%%%note there will be 0/0 = NaN, but the SZA criteria will skip these.
XJ1 = J11./J12;
XJ2 = J21./J22;
lnXJ1 = log(XJ1);
lnXJ2 = log(XJ2);

for i=1:NX
 for j=1:NY
  for b=1:Blocks
   if block(i,j,b) == 1
    for t=1:24
       ishift = floor(t*NX/24.);
       ii = i + ishift;
     if ii>NX
       ii = ii-NX;
     end
     if cosy(ii,j) > SZA0
      for l=1:NP
        k = max(1,min(9,floor(plevel(i,j,l)/100)));
       Jk1 = min(201, max(1, 101+floor(100*lnXJ1(i,j,l,t)+0.5)));
       Jk2 = min(201, max(1, 101+floor(100*lnXJ2(i,j,l,t)+0.5)));
       pdMJ(Jk1,k,1,b,model) = pdMJ(Jk1,k,1,b,model) + 1;
       pdMJ(Jk2,k,2,b,model) = pdMJ(Jk2,k,2,b,model) + 1;
      end
     end
    end
   end
  end
 end
end

for b=1:Blocks
 for j=1:2
  for k=1:9
    sumMJ = sum(pdMJ(:,k,j,b,model));
   if sumMJ > 0
    pdMJ(:,k,j,b,model) = pdMJ(:,k,j,b,model)/sumMJ;
   end
    pdMJbot(k,j,b,model) = sum(pdMJ(1:50,k,j,b,model));
    pdMJtop(k,j,b,model) = sum(pdMJ(152:201,k,j,b,model));
  end
 end
end

for k3=1:3
 for j=1:2
  for b=1:Blocks
    sumJ(1:201) = 0;
    cntJ = 0;
   for k = Plev1(k3):Plev2(k3)
     sumMJ = sum(pdMJ(:,k,j,b,model));
   if sumMJ > 0
    for n=1:201
     sumJ(n) = sumJ(n) + pdMJ(n,k,j,b,model);
    end
     cntJ = cntJ + 1;
   end
   end
    for n=1:201
     JXstat(n,k3,j,b,model) = sumJ(n)/cntJ;
    end
  end
 end
end


%%%%% end of loop over all models
end
%%%%%%