close all;
clear all;

load conus;
S = shaperead('landareas','UseGeoCoords',true);

%do the analysis over boxrad * boxrad area using data within (boxrad+1)x(boxrad+1)area

load('inputs.dat');
xc = inputs(1);
yc = inputs(2);
boxrad_lon = inputs(3);
boxrad_lat = inputs(4);
sigo = inputs(5);
sigb = inputs(6);
rad = inputs(7);
radz = inputs(8);
subNo = inputs(9);
radius = 6370;
rad2 = rad * rad;
radz2 = radz * radz;

% set obs threshold for background
thrhdb = 10;
% set obs threshold for anlysis
thrhda = 6;

% set plotting range
xmax = xc + boxrad_lon;
xmin = xc - boxrad_lon;
ymax = yc + boxrad_lat;
ymin = yc - boxrad_lat;

xmaxb = xmax;
xminb = xmin;
ymaxb = ymax;
yminb = ymin;

%read background file temperature
filename_temp = strcat('./bkg_data.dat');
load (filename_temp)

% read data into individual arrays
numxb=0;
numbkgpts = numel(bkg_data(:,1));
for k=1:numbkgpts;
    if ( (bkg_data(k,2) >= xminb) && (bkg_data(k,2) <= xmaxb) && (bkg_data(k,1) >= yminb) && (bkg_data(k,1) <= ymaxb))
        numxb = numxb+1;
        xb_lat(numxb,1) = bkg_data(k,1);
        xb_lon(numxb,1) = bkg_data(k,2);
        xb_elev(numxb,1) = bkg_data(k,3);
%         xbt(numxb,1) = bkg_data(k,4);
        xbt(numxb,1) = 0;
    end
end

clear bkg_data input numxb;

numpts = size(xb_lat,1);
pb=zeros(numpts,numpts);
for j=1:numpts
    dx = radius .* cos(pi .*xb_lat ./180.) .* pi .* (xb_lon - xb_lon(j))./ 180.;
    dy = radius .* pi .* (xb_lat - xb_lat(j)) ./ 180.;
    dz = xb_elev - xb_elev(j); 
    r2 = dx .* dx + dy .* dy;
    z2 =  (dz .* dz);
    pb(:,j) = sigb .* (exp(-r2./rad2).*exp(-z2/radz2));
end
clear dx dy dz r2 z2;

% load observation data
filename_obs = strcat('./obs.dat');
[o_yr,o_mn,o_dy,o_hr,o_min,o_id,o_mnet,o_lat,o_lon,o_elev,o_t,o_r,o_s,o_d] = textread(filename_obs,'%f %f %f %f %f %s %f %f %f %f %f %f %f %f');

% find the size of the array
numobst = size(o_lat,1);


o_t=ones(numobst,1);

jj=0;
ii=0;
for j=1:numobst;
    if(o_t(j)== -9999)  %remove missing data
    else
        % compute min_val and min_index
        yo_lat=o_lat(j);
        yo_lon=o_lon(j);
        dx = radius .* cos(pi .* yo_lat ./ 180.) .* pi .*(yo_lon-xb_lon)./ 180.;
        dy = radius .* pi .* (yo_lat - xb_lat) ./ 180.;
        dist = sqrt(dx.*dx + dy.*dy);
        [min_val,min_index] = min(dist);
        % do simple check against background to exclude potentially bad obs
%         dfob=5*(o_t(j)-32)/9-xbt(min_index);
        dfob=o_t(j)-xbt(min_index);
        if(abs(dfob)<thrhdb)
            jj=jj+1;
            index(jj) = min_index;
            h(jj,:)=zeros(1,numpts);
            h(jj,min_index) = 1;
            ob_yr(jj)=o_yr(j);
            ob_mn(jj)=o_mn(j);
            ob_dy(jj)=o_dy(j);
            ob_hr(jj)=o_hr(j);
            ob_min(jj)=o_min(j);
            ob_id(jj)=o_id(j);
            ob_mnet(jj)=o_mnet(j);
            ob_lat(jj)=o_lat(j);
            ob_lon(jj)=o_lon(j);
            ob_elev(jj)=o_elev(j);
%             ob_t(jj)=5*(o_t(j)-32)/9;
            ob_t(jj)=o_t(j);

            %text coordinates
            latb(jj,ob_mnet(jj))=ob_lat(jj);
            lonb(jj,ob_mnet(jj))=ob_lon(jj);
            drb(jj,1)=ob_t(jj)-xbt(min_index);
            dmtb(jj,ob_mnet(jj))=ob_t(jj)-xbt(min_index);
            tmb(jj,ob_mnet(jj))=ob_t(jj);
            
            if((o_lat(j)>=ymin)&&(o_lat(j)<=ymax)&&(o_lon(j)>=xmin)&&(o_lon(j)<=xmax));
                % these are obs in data denial category
                 ii=ii+1;
                 oo_bin(ii)=jj;
            end
        end
    end
end

numobs = jj;
numobsi = ii;
numob=numobs;
% 3Dvar conjugate gradient solutions
% temp=pb*h'
temp = zeros(numpts,numob);
for i=1:numpts;
    for j=1:numob;
        temp(i,j)=pb(i,index(j));
    end
end
% the linear system is equivalent to the following: 
pbt = zeros(numpts,numpts);
pbhpo = zeros(numpts,1);
pbt = sigo * eye(numpts,numpts) + temp * h;
pbhpo = temp * drb;
x3dt = xbt + gmres(pbt,pbhpo,1000,0.001);
clear temp pbt pbhpo;
    
% save control analysis
x3dtc = x3dt;

filename = strcat('scaler_',num2str(yc),'_',num2str(xc),'.dat');
fid=fopen(char(filename),'w');
for i=1:numpts 
    fprintf(fid,'%f %f %f %f',xb_lat(i), xb_lon(i), xb_elev(i),x3dtc(i));
    fprintf(fid,'\n');
end

% move the obs that does not meet threshold for analysis

mm=0;
for j=1:numobst;
    if(o_t(j)== -9999)  %remove missing data
    else
        % compute min_val and min_index
        yo_lat=o_lat(j);
        yo_lon=o_lon(j);
        dx = radius .* cos(pi .* yo_lat ./ 180.) .* pi .*(yo_lon-xb_lon)./ 180.;
        dy = radius .* pi .* (yo_lat - xb_lat) ./ 180.;
        dist = sqrt(dx.*dx + dy.*dy);
        [min_val,min_index] = min(dist);
        % do simple check against analysis to exclude potentially bad obs
%         dfob=5*(o_t(j)-32)/9-x3dtc(min_index);
        dfob=o_t(j)-x3dtc(min_index);
        if(abs(dfob)<thrhda)
            mm=mm+1;
            oa_yr(mm)=o_yr(j);
            oa_mn(mm)=o_mn(j);
            oa_dy(mm)=o_dy(j);
            oa_hr(mm)=o_hr(j);
            oa_min(mm)=o_min(j);
            oa_id(mm)=o_id(j);
            oa_mnet(mm)=o_mnet(j);
            oa_lat(mm)=o_lat(j);
            oa_lon(mm)=o_lon(j);
            oa_elev(mm)=o_elev(j);
%             oa_t(mm)=5*(o_t(j)-32)/9;
            oa_t(mm)=o_t(j);


            %text coordinates
            lata(mm,oa_mnet(mm))=oa_lat(mm);
            lona(mm,oa_mnet(mm))=oa_lon(mm);
            dra(mm,1)=oa_t(mm)-x3dtc(min_index);
            dmta(mm,oa_mnet(mm))=oa_t(mm)-x3dtc(min_index);
            tma(mm,oa_mnet(mm))=oa_t(mm);
        end
    end
end

clear obs o_yr o_mn o_dy o_hr o_min o_id o_mnet o_lat o_lon o_elev o_t o_r o_s o_d o_u o_v;


outfigdir = './';          % figure output directory, needs / on end
cmap = [0.2235,0.1137,0.5961;0.1529,0.2118,0.6235;0.1059,0.3294,0.7216;0.0588,0.4667,0.7843;0.0353,0.6667,0.8549;0.0078,0.7529,0.9412; ...
        0,0.8118,0.9176;0.0039,0.8627,0.7725;0.0039,0.8549,0.5647;0,0.8039,0.2706;0.0275,0.7529,0.0863;0.1137,0.7608,0.0039; ...
       0.302,0.7922,0.0078; 0.5647,0.8235,0.0118;0.8235,0.8745,0.0196;0.9412,0.9216,0.0118;0.9922,0.9843,0.0039;1,0.9216,0; ...
       1,0.8235,0;0.9961,0.702,0; 0.9882,0.5059,0;0.9647,0.3098,0;0.9333,0.2235,0;0.8824,0.1843,0.0078;0.8353,0.1451,0.0078; ...
       0.7843,0.1059,0.0078;0.7294,0.0745,0.0078;0.6824,0.0549,0.0078;0.6353,0.0353,0.0078;0.5098,0.0353,0.0078;];
% orange/purple map for increments
% 20 colors
opcmap = [0.3529,0.05490,0.4745;0.3765,0.06670,0.4980;0.4353,0.09410,0.5529;0.5098,0.1529,0.6157;0.5922,0.2353,0.6863;0.6784,0.3412,0.7569;...
          0.7765,0.4941,0.8314;0.8627,0.6667,0.8980;0.9529,0.8784,0.9647;1,1,1;1,1,1;0.9961,0.9765,0.7922;0.9922,0.9294,0.4824;0.9882,0.8863,0.2863;...
          0.9843,0.8392,0.1569;0.9804,0.7843,0.07450;0.9765,0.7294,0.03530;0.9725,0.6784,0.01570;0.9686,0.6353,0.007800;0.9647,0.6157,0.0039;];

xi = xmin:0.05:xmax;
yi = ymin:0.05:ymax;
xib = xminb:0.04:xmaxb;
yib = yminb:0.04:ymaxb;
[XI,YI] = meshgrid(xi,yi);
[XIB,YIB] = meshgrid(xib,yib);

% plot background temperature over sub domain
ZZIB = griddata(xb_lon,xb_lat,xbt,XIB,YIB,'cubic');
clear numobs h v1 v2 t scl;
numobs=size(latb,1);
numoas=size(lata,1);

% %--------------------------------------------------------------------------
% % plot temperature increments over subdomain
% ZZII = griddata(xb_lon,xb_lat,x3dtc-xbt,XI,YI);
% v1=0;
% v2=1;
% v = v1:0.1:v2;
% h=0.1;
% t=(v2-v1)/18;
% scl = v1:t:v2;
% 
% figure(2);
% box off;
% axis off;
% axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'PLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'MLabelRound',-1);
% contourfm(YI,XI,ZZII,v);
% shading flat;
% colormap(opcmap);
% cbh=colorbar;
% caxis([v1 v2]);
% set(cbh,'YTick',v1:h:v2);
% [C,h]=contourm(YIB,XIB,ZI,[min_elev:H:max_elev],'k--');
% % clabel(C,h,[500:500:3000],'FontName','Times-Roman','FontSize',8,'Color','k','Rotation',0,'FontWeight','Bold');
% % textm(latb(:,1),lonb(:,1),num2str(dmtb(:,1),'%3.0f'),'FontName','Times-Roman','FontSize',10,'Color','b','HorizontalAlignment','center','Clipping','on');
% % textm(latb(:,2),lonb(:,2),num2str(dmtb(:,2),'%3.0f'),'FontName','Times-Roman','FontSize',10,'Color','g','HorizontalAlignment','center','Clipping','on');
% title('2D-Var Innovations - Temperature');
% 
% for k=1:numobs;
%     plotm(latb(k,1),lonb(k,1),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     % determine the marker color based on the temperature scale
%     lngthscl = length(scl);
%     % NWS    
%     if ( dmtb(k,1) <= scl(1) )
%         plotm(latb(k,1),lonb(k,1),'o','MarkerEdgeColor',opcmap(1,:),'MarkerFaceColor',opcmap(1,:),'MarkerSize',5);
%     elseif ( dmtb(k,1) >= scl(lngthscl) )
%         plotm(latb(k,1),lonb(k,1),'o','MarkerEdgeColor',opcmap(lngthscl,:),'MarkerFaceColor',opcmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (dmtb(k,1) > scl(i-1)) && dmtb(k,1) <= scl(i))
%                 plotm(latb(k,1),lonb(k,1),'o','MarkerEdgeColor',opcmap(i,:),'MarkerFaceColor',opcmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
%     % RAWS
%     if ( dmtb(k,2) <= scl(1) )
%         plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(1,:),'MarkerFaceColor',opcmap(1,:),'MarkerSize',5);
%     elseif ( dmtb(k,2) >= scl(lngthscl) )
%         plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(lngthscl,:),'MarkerFaceColor',opcmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (dmtb(k,2) > scl(i-1)) && dmtb(k,2) <= scl(i))
%                 plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(i,:),'MarkerFaceColor',opcmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
% end;
% 
% if (subNo==1||subNo==11||subNo==19||subNo==20||subNo==32||subNo==33||subNo==31||subNo==42||subNo==49||subNo==46||subNo==47||subNo==48)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==2||subNo==3||subNo==7||subNo==8||subNo==9||subNo==10||subNo==18)
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==4||subNo==5||subNo==44)
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==45)
%     disp('No lines');
% elseif(subNo==50||subNo==51)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
% else    
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% end
% 
% filename = strcat(outfigdir,'scaler_inct.eps');
% saveas(gcf,filename,'epsc2');
% 
% %--------------------------------------------------------------------------
% % set scaler for analysis 
% v1=0;
% v2=1;
% v=v1:0.1:v2;
% h=0.1;
% t=(v2-v1)/25;
% scl = v1:t:v2;
% ZZV = griddata(xb_lon,xb_lat,x3dtc,XI,YI,'cubic');
% 
% figure(3)
% box off;
% axis off;
% axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'MLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'PLabelRound',-1);
% contourfm(YI,XI,ZZV,v);
% shading flat;
% colormap(cmap);
% cbh=colorbar;
% caxis([v1 v2]);
% set(cbh,'YTick',v1:h:v2);
% [C,h]=contourm(YIB,XIB,ZI,[min_elev:H:max_elev],'k--');
% % clabel(C,h,[500:500:3000],'FontName','Times-Roman','FontSize',8,'Color','k','Rotation',0,'FontWeight','Bold');
% % textm(lata(:,1),lona(:,1),num2str(tma(:,1),'%3.0f'),'FontName','Times-Roman','FontSize',10,'Color','b','HorizontalAlignment','center','Clipping','on');
% % textm(lata(:,2),lona(:,2),num2str(tma(:,2),'%3.0f'),'FontName','Times-Roman','FontSize',10,'Color','g','HorizontalAlignment','center','Clipping','on');
% title('2D-Var Analysis  - Temperature');
% for k=1:numoas;
%     plotm(lata(k,1),lona(k,1),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     % determine the marker color based on the temperature scale
%     lngthscl = length(scl);
%     % NWS
%     if ( tma(k,1) <= scl(1) )
%         plotm(lata(k,1),lona(k,1),'o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',5);
%     elseif ( tma(k,1) >= scl(lngthscl) )
%         plotm(lata(k,1),lona(k,1),'o','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (tma(k,1) > scl(i-1)) && tma(k,1) <= scl(i))
%                 plotm(lata(k,1),lona(k,1),'o','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
%     %RAWS
%     if ( tma(k,2) <= scl(1) )
%         plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',5);
%     elseif ( tma(k,2) >= scl(lngthscl) )
%         plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (tma(k,2) > scl(i-1)) && tma(k,2) <= scl(i))
%                 plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
% end;
% 
% if (subNo==1||subNo==11||subNo==19||subNo==20||subNo==32||subNo==33||subNo==31||subNo==42||subNo==49||subNo==46||subNo==47||subNo==48)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==2||subNo==3||subNo==7||subNo==8||subNo==9||subNo==10||subNo==18)
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==4||subNo==5||subNo==44)
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==45)
%     disp('No lines');
% elseif(subNo==50||subNo==51)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
% else    
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% end
% 
% filename = strcat(outfigdir,'scaler_anlt.eps');
% saveas(gcf,filename,'epsc2');


%--------------------------------------------------------------------------
% % plot temperature increments over subdomain
% opcmap=zeros(25,3);
% for i=1:26
%     opcmap(i,:)=(26-i)*0.04;
% end

elev_scl=[50,500,1500,2500];
ZI = griddata(xb_lon,xb_lat,xb_elev,XIB,YIB,'cubic');

ZZII = griddata(xb_lon,xb_lat,x3dtc-xbt,XI,YI);
% 
% v1=0;
% v2=1;
% h=0.1;
% t=(v2-v1)/24;
% v = v1:t:v2;
% scl = v1:t:v2;
% 
% figure(2);
% box off;
% axis off;
% axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'PLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'MLabelRound',-1);
% contourfm(YI,XI,ZZII,v);
% shading flat;
% colormap(opcmap);
% cbh=colorbar;
% caxis([v1 v2]);
% set(cbh,'YTick',v1:h:v2);
% [C,h]=contourm(YIB,XIB,ZI,elev_scl,'k--');
% clabel(C,h,elev_scl,'FontName','Times-Roman','FontSize',6,'Color','k','Rotation',0);
% 
% for k=1:numobs;
%     plotm(latb(k,1),lonb(k,1),'s','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',5);
%     plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     % determine the marker color based on the temperature scale
%     lngthscl = length(scl);
%     % NWS    
%     if ( dmtb(k,1) <= scl(1) )
%         plotm(latb(k,1),lonb(k,1),'s','MarkerEdgeColor',opcmap(1,:),'MarkerFaceColor',opcmap(1,:),'MarkerSize',3.5);
%     elseif ( dmtb(k,1) >= scl(lngthscl) )
%         plotm(latb(k,1),lonb(k,1),'s','MarkerEdgeColor',opcmap(lngthscl,:),'MarkerFaceColor',opcmap(lngthscl,:),'MarkerSize',3.5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (dmtb(k,1) > scl(i-1)) && dmtb(k,1) <= scl(i))
%                 plotm(latb(k,1),lonb(k,1),'s','MarkerEdgeColor',opcmap(i,:),'MarkerFaceColor',opcmap(i,:),'MarkerSize',3.5);
%             end 
%         end;
%     end
%     % RAWS
%     if ( dmtb(k,2) <= scl(1) )
%         plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(1,:),'MarkerFaceColor',opcmap(1,:),'MarkerSize',5);
%     elseif ( dmtb(k,2) >= scl(lngthscl) )
%         plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(lngthscl,:),'MarkerFaceColor',opcmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (dmtb(k,2) > scl(i-1)) && dmtb(k,2) <= scl(i))
%                 plotm(latb(k,2),lonb(k,2),'o','MarkerEdgeColor',opcmap(i,:),'MarkerFaceColor',opcmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
% end;
% 
% if (subNo==1||subNo==11||subNo==19||subNo==20||subNo==32||subNo==33||subNo==31||subNo==42||subNo==49||subNo==46||subNo==47||subNo==48)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==2||subNo==3||subNo==7||subNo==8||subNo==9||subNo==10||subNo==18)
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==4||subNo==5||subNo==44)
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==45)
%     disp('No lines');
% elseif(subNo==50||subNo==51)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
% else    
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% end
% 
% filename = strcat(outfigdir,'scaler_inct_gray.eps');
% saveas(gcf,filename,'epsc2');

%--------------------------------------------------------------------------
% set scaler for analysis 
cmap=zeros(25,3);
for i=1:26
    cmap(i,:)=1-(i-1)*0.03;
end


v1=0;
v2=1;
h=0.1;
t=(v2-v1)/24;
v=v1:t:v2;
scl = v1:t:v2;
ZZV = griddata(xb_lon,xb_lat,x3dtc,XI,YI,'cubic');

% figure(3)
% box off;
% axis off;
% axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'MLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'PLabelRound',-1);
% contourfm(YI,XI,ZZV,v);
% shading flat;
% colormap(cmap);
% cbh=colorbar;
% caxis([v1 v2]);
% set(cbh,'YTick',v1:h:v2);
% [C,h]=contourm(YIB,XIB,ZI,elev_scl,'k--');
% clabel(C,h,elev_scl,'FontName','Times-Roman','FontSize',6,'Color','k','Rotation',0);
% 
% for k=1:numoas;
%     plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',6);
%     plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
%     % determine the marker color based on the temperature scale
%     lngthscl = length(scl);
%     % NWS
%     if ( tma(k,1) <= scl(1) )
%         plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',3.5);
%     elseif ( tma(k,1) >= scl(lngthscl) )
%         plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',3.5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (tma(k,1) > scl(i-1)) && tma(k,1) <= scl(i))
%                 plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',3.5);
%             end 
%         end;
%     end
%     %RAWS
%     if ( tma(k,2) <= scl(1) )
%         plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',5);
%     elseif ( tma(k,2) >= scl(lngthscl) )
%         plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',5);
%     else
%         % loop through the colors in the scale
%         for i=2:lngthscl;
%             if ( (tma(k,2) > scl(i-1)) && tma(k,2) <= scl(i))
%                 plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',5);
%             end 
%         end;
%     end
% end;
% 
% if (subNo==1||subNo==11||subNo==19||subNo==20||subNo==32||subNo==33||subNo==31||subNo==42||subNo==49||subNo==46||subNo==47||subNo==48)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==2||subNo==3||subNo==7||subNo==8||subNo==9||subNo==10||subNo==18)
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==4||subNo==5||subNo==44)
%     geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
% elseif(subNo==45)
%     disp('No lines');
% elseif(subNo==50||subNo==51)
%     geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
% else    
%     geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
% end
% 
% filename = strcat(outfigdir,'scaler_anlt_gray.eps');
% saveas(gcf,filename,'epsc2');


v1=0;
v2=3400;
t=(v2-v1)/24;
v=v1:t:v2;
figure(31)
box off;
axis off;
axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'MLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'PLabelRound',-1);
contourfm(YIB,XIB,ZI,v);
shading flat;
colormap(cmap);
cbh=colorbar;
caxis([0 3000]);
set(cbh,'YTick',v1:500:v2);
[C,h]=contourm(YIB,XIB,ZI,elev_scl,'k--');
clabel(C,h,elev_scl,'FontName','Times-Roman','FontSize',6,'Color','k','Rotation',0);
for k=1:numoas;
    plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',6);
    plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',7);
    % determine the marker color based on the temperature scale
    lngthscl = length(scl);
    % NWS
    if ( tma(k,1) <= scl(1) )
        plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',3.5);
    elseif ( tma(k,1) >= scl(lngthscl) )
        plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',3.5);
    else
        % loop through the colors in the scale
        for i=2:lngthscl;
            if ( (tma(k,1) > scl(i-1)) && tma(k,1) <= scl(i))
                plotm(lata(k,1),lona(k,1),'s','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',3.5);
            end 
        end;
    end
    %RAWS
    if ( tma(k,2) <= scl(1) )
        plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',cmap(1,:),'MarkerSize',5);
    elseif ( tma(k,2) >= scl(lngthscl) )
        plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(lngthscl,:),'MarkerFaceColor',cmap(lngthscl,:),'MarkerSize',5);
    else
        % loop through the colors in the scale
        for i=2:lngthscl;
            if ( (tma(k,2) > scl(i-1)) && tma(k,2) <= scl(i))
                plotm(lata(k,2),lona(k,2),'o','MarkerEdgeColor',cmap(i,:),'MarkerFaceColor',cmap(i,:),'MarkerSize',5);
            end 
        end;
    end
    
end;
plotm(39.75972,-121.6208,'p','MarkerEdgeColor','k','MarkerFaceColor','w','MarkerSize',11);
textm(39.32,-120.68,'KBLU','color','w','FontSize',9);
if (subNo==1||subNo==11||subNo==19||subNo==20||subNo==32||subNo==33||subNo==31||subNo==42||subNo==49||subNo==46||subNo==47||subNo==48)
    geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
    geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
elseif(subNo==2||subNo==3||subNo==7||subNo==8||subNo==9||subNo==10||subNo==18)
    geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
    geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
elseif(subNo==4||subNo==5||subNo==44)
    geoshow(uslat,uslon,'DisplayType','line','Color','k','LineWidth',1);
elseif(subNo==45)
    disp('No lines');
elseif(subNo==50||subNo==51)
    geoshow([S.Lat], [S.Lon],'Color','k','LineWidth',1);
else    
    geoshow(statelat,statelon,'DisplayType','line','Color','k','LineWidth',1);
end

filename = strcat(outfigdir,'scaler_terrain_gray1.eps');
saveas(gcf,filename,'epsc2');