%% changed from tig_cmp2_plt8.m
% Upload with publications.
    close all; clear;
    dash_dir = '/glade/work/ldong/DASH/Lin24_CIMI_HF/';
    ncfn = 'CIMI_HF_fig2.nc';
    ncname = [dash_dir,ncfn];
    disp(['Loading data from ',ncname]);
    
    ncid = netcdf.open(ncname,'NC_NOWRITE');
    info = ncinfo(ncname);
    varnames = {}; % store all available variable names.
    for ii=1:size(info.Variables,2)
        varnames{ii} = info.Variables(ii).Name;
        disp([num2str(ii,'%02d'),' ',info.Variables(ii).Name]);
    end

    for ii=1:length(varnames)
        varid = netcdf.inqVarID(ncid,varnames{ii});
        eval([varnames{ii},' = netcdf.getVar(ncid,varid);']);
    end
    
    %%

        [~,i1] = max((dCER(:))); 
        [ilonm,ilatm] = ind2sub(size(dCER),i1);
        
        dCER_plt  = dCER(isort,:);    % plot against LON/LAT
        [~,i1] = max((dCER_plt(:))); 
        [ilon,ilat] = ind2sub(size(dCER_plt),i1);
        % [LON(ilonm,ilatm),LAT(ilonm,ilatm),LON(isort(ilon),ilatm),LAT(isort(ilon),ilatm)]
        % % to verify the lat/lon for max CER is corret in the sorted
        % array.
%         [~,ilonp] = min(abs(mod(LON(:,1)/15+24,24)-LT(ilon,ilat)));
%         [ilonm,ilatm,LT(ilonm,ilatm),LON(ilonm,ilatm),LAT(ilonm,ilatm),LON(ilonp,1)/15,ilonp]

        % alt_plt = squeeze(ZG1(isort(ilon),ilat,:))*1e-5;
        varplt1d1s = {squeeze((Te1(isort(ilon),ilat,:))), squeeze(log10(Ne1(isort(ilon),ilat,:))), ...
            squeeze(VER1(isort(ilon),ilat,:))};
        varplt1d2s = {squeeze((Te2(isort(ilon),ilat,:))), squeeze(log10(Ne2(isort(ilon),ilat,:))), ...
            squeeze(VER2(isort(ilon),ilat,:))};
        varplt2ds = {dTe(isort,:), dNe(isort,:), dCER(isort,:)};

        dfelim = [8,10];
        dtelim = [-1,1]*1000;
        dnelim = [-1,1]*1e4;
        dcerlim = [-1,1]*200;
        dTNlim = [-1,1]*1e0;
        clims = {dtelim, dnelim, dcerlim, dTNlim};
        dtemap= '*RdBu';
        cmaps = {dtemap,dtemap,dtemap,dtemap};
        cttls = {'\Delta T_e [K]','\Delta N_e [/cc]','\Delta CER [R]','\Delta T_N [K]'};

        Telim = [1000,5000];
        Nelim = [4,4.8]+0.2;
        verlim= [0,32];
        TNlim = [0,1000];
        Zlim = [90,600];
        x1dlims = {Telim,Nelim,verlim,TNlim};
        cl1 = 'k';
        cl2 = 'r';
        xlabs = {'T_e [K]','N_e log_{10}[/cc]','VER [Photons/cm^3/s]','T_N [K]'};
        
        % ttl = {['Pressure level ',num2str(lev(ialt_smp),'%.1f'),' (~',num2str(mean(oneD(ZG1(:,:,ialt_smp)*1e-5)),'%.0f'),' km)']};
        ttl = 'Pressure level 2.1 (~293 km)';
        
        slt2lon = 15;
        latlim = [-1,1]*90;
        sltlim = [0,24];
        lonlim = sltlim*slt2lon-180;
        lattick = -80:20:80;
        slttick = 2:2:22;
        lontick = slttick*slt2lon-180;
        sltticklab = {};
        for ii=1:length(slttick)
            sltticklab{ii} = num2str(slttick(ii),'%02d');
        end
        sltticklab;
        cl_grds = {'k','k','k','k'};
        cl_llb = [0.4660    0.6740    0.1880]; %cl_gre;
        
        do_cont = 0;
        panlabs = {'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)'};
        cl_pans = {'k','k','k','k','k','k','k','k'};
        xt1 = 188;
        yt1 = 135;
        xpan1 = -177;
        ypan1 = 70;
        xpanr2 = 0.02;
        ypanr2 = 0.95;
        strdate = '2013-06-01T01-00-00';
        %
        
        hps = {}; cls = {};
        font = 14; lwid=2;
        pngname = [dash_dir,'CIMI_HF_fig2'];
        disp(['Plotting ',pngname]);
        nrow = 3; ncol = 2; fpos = [1,39,1200,900];
        fig = figure('position',fpos,'visible','off');
    for irow=1:nrow
        isub = (irow-1)*ncol+1;
        hps{isub} = subplot(nrow,ncol,isub);
        hold on; box on; grid on; axis equal tight;
        set(gca,'Fontsize',font,'Linewidth',lwid);
        colormap(hps{isub},brewermap([],cmaps{irow}));
        pcolor(LON,LAT,varplt2ds{irow}); shading flat; 
        cls{isub} = colorbar; 
        if isub>0
            plot(LON(ilon,ilat),LAT(ilon,ilat),'m*','Linewidth',lwid*1.5,'color','m','Markersize',10)
        end
        caxis(clims{irow});
        if do_cont
            contour(LON,LAT,varplt2ds{irow},varlvls{irow},'m','ShowText','on'); 
        end
        plot(aurllbN(:,1),aurllbN(:,2),'Linewidth',lwid,'color',cl_llb);
        plot(aurllbS(:,1),aurllbS(:,2),'Linewidth',lwid,'color',cl_llb);
        % % contour(LON,LAT,geng1(isort,:),[0.1,1.0],'c'); 
        % contour(LON,LAT,gflx1(isort,:),[1,1]*nflx0,'g'); 
        xlim(lonlim);
        ylim(latlim);
        if irow==nrow
            xlabel('SLT [Hour]');
        end
        ylabel('Geo Lat [deg]');
        set(gca,'xtick',lontick,'ytick',lattick,'xticklabel',sltticklab);
        plot_grid(lontick,lattick,cl_grds{irow});
        if isub==1
            title(ttl);
            htt = text(xt1,yt1,strdate,'Fontsize',font+1,'Fontweight','Bold');
        end
        ylabel(cls{isub},cttls{irow});
        text(xpan1,ypan1,panlabs{isub},'Fontsize',font,'color',cl_pans{isub});
    end
        
    for irow=1:nrow
        isub = (irow-1)*ncol+2;
        hps{isub} = subplot(nrow,ncol,isub);
        hold on; box on; grid on; axis tight;
        set(gca,'Fontsize',font,'Linewidth',lwid);
        plot(varplt1d1s{irow},alt_plt,'color',cl1,'Linewidth',lwid);
        plot(varplt1d2s{irow},alt_plt,'color',cl2,'Linewidth',lwid);
        ylim(Zlim);
        if irow<nrow
            xlim(x1dlims{irow});
        end
        xlabel(xlabs{irow}); 
        ylabel('Alt [km]');
        if irow==1
            legend({'Baseline','CIMI Heat'},'Location','SouthEast');
            title({['Lon=',num2str(LON(isort(ilon),ilat),'%.1f'),'^\circ, Lat=',num2str(LAT(isort(ilon),ilat),'%.1f'),'^\circ']});
        end
        xlim1 = get(gca,'xlim');
        ylim1 = get(gca,'ylim');
        text(xlim1(1)+xpanr2*(xlim1(2)-xlim1(1)),ylim1(1)+ypanr2*(ylim1(2)-ylim1(1)),panlabs{isub},'Fontsize',font);
    end

        % for isub=4:5
        %     hpos = get(hps{isub},'position');
        %     set(hps{isub},'position',hpos+[0.02,0,0,0]);
        % end
        
        % set(gcf,'renderer','painters');
        print(fig,pngname,'-dpng');
        close(fig);
        
        %%
        save_nc = 0;
        
        if save_nc
            %% DMSP_tplt0,log10(DMSP_chen0),log10(DMSP_eden0)
            disp('Do you really want to overwrite the nc files?');
            pause
            dash_dir = '/glade/work/ldong/DASH/Lin24_CIMI_HF/';
            if(~exist(dash_dir,'dir')); mkdir(dash_dir); end
    
            ncfn = 'CIMI_HF_fig2.nc';
            ncname = [dash_dir,ncfn];
            disp(['Writing data to ',ncname]);
    
            % Open a netcdf file:
            ncid = netcdf.create(ncname,'NC_WRITE');
            varnames = {'LON','LAT','dTe','dNe','dCER','isort','ilon','ilat','aurllbN','aurllbS',...
                'Te1','Ne1','VER1','Te2','Ne2','VER2','alt_plt'};
            run("save_nc_file_cell.m");
        end