% um_plot.m
% Basic plotting program on model domain
% gnp sorted out the time problem.
% Creates a 2-d contour plot (using contour, contourf or pcolor)
% one or two variables plus vectors can be plotted
% Cross-sections chosen via scntype: 1=x-y, 2=x-z, 3=y-z , 4=interpolated
% For x-y plots the aspect ratio is adjusted so dx=dy and the
% coastline (via lsm) and a lat/lon grid are plotted
% use after um_read_nc; 
% calls um_set_field and um_add_vectors and um_time_test and um_level_test
% written by Ian Renfrew, Aug 2006
% modified by Gudrun Nina Petersen and Stephen Outten, 2007-2009
% based on code provide by Jeff Chagnon
%------------------------------------------------
close all; clear field1 plotfield conlevs conlevs_labels slice
%-----------------------------------------------------------------
% set up main plotting parameters here
%------------------------------------------------------------------
iplot = 2;          % 1=contour, 2=contourf, 3=pcolor, 4=contour+con, 5=contourf+con, 6=pcolor+con 
scntype = 4;        % section type: 1 = x-y, 2 = y-z, 3 = x-z, 4 = interpolated slice
iheight = 1;        % defines height for y-z, x-z plots: 0=hybrid height, 1=geometric height
height_limit = 4;  % top of plotting domain for x-z or y-z plots (km)
ilev1 = 1;          % set vertical level for x-y plot; for single levels ilev is set to 1 later; level 12 for OSTIA_ice+sst best level
ilev2 = ilev1;      % Set as different from ilev1 if necessary
iy = 55;            % set model grid index for x-z plots
ix = 164;           % set model grid index for y-z plots
itime = [1:1];        % set time index: scalar or array e.g. [1:1:5]
%%itime2 = 0;         % 0=one var or time array same for both variables. 
iprint = 0;         % print to a file if iprint =1 (ps) or 2 (jpg) or 3 (psc2)
icolormap = 0;      % colormap: 0=jet, 1= yellow-red, 2= white-black, 
icolorbar = 0;      % plot colorbar if 1
ivec = 0;           % plot vectors if ivec~=0; plot every nth in i-dir 
jvec = 8;           % plot every nth in j-direction
AX=71; AY=57; BX=62; BY=58;  % ************ defines point A and B for interpolated slice e.g.(RTJ - DS1-4) itime = 2
slice_limit = 100;  % defines length of interpolated slice RTJ slice1=100, slice2=269, slice3=119, s1=430, s2=595, s3=437, s4=403
dl=1;     % ******** sets resolution of interpolated slice e.g. 1km 
zoom_plot =0;      % turns the zoom on and off (1 or 0)
zooma =1; zoomb =100; zoomc=1; zoomd=100; % defines the two points used for the zoom (zooma,zoomc) and (zoomb,zoomd)
plot_over_orog = 1;    % 0=variable not ploted over land, 1=variable plotted over land
plot_vect_over_orog = 1;    % 0=vectors not ploted over land, 1=vectors plotted over land
%-----------------------------------------------------------------
% vertical coordinate prescribed via vtype =
%    1= hheight_rho; 2=hheight_theta; 3=theta_levels; 4=p_levels; 5=single
%
% Set time in hours, set as t, t_1, t_2, t_3, etc as appropriate
%
% Variable names are defined in um_read_nc:
% atmospheric: model levels, t, T6Hourly:    1=u, 2=v, 3=theta, 4=q, 5=w, 6=pv
% atmospheric: model levels, t_3, T6Hourly:  7=qcl, 8=qcf
% atmospheric: theta levels, t, T6Hourly:    9=pv_theta,
% atmospheric: pressure levels, t, T6Hourly: 10=gph
% single level: t_1, T0 (time=0 only):       11=lsm, 12=orography
% single level: t, T6Hourly:                 13=temp_sfc, 14=mslp
% single level: t_3, T6Hourly: 15=temp_1pt5m, 16=q_1pt5m, 17=u_10m, 18=v_10m, 19=wind_speed_10m
% single level: t_3, T6Hourly: 20=low_cloud, 21=med_cloud, 22=high_cloud
% single level: t_2, T6Mean6: 
%  23=shf_sfc, 24=lhf_sfc, 25=taux_sfc, 26=tauy_sfc, 27=tau_sfc
%  28=lsrain, 29=lssnow, 30=cvrain, 31=cvsnow, 32=wind_speed, 33=p_levels1,
%  34=wind_dir, 35=ice_conc
% Add new diagnostics here......
%-----------------------------------------------------------------
ivar = 3;             % SET VARIABLES HERE; contour levels can be adjusted in um_set_field
ivar2 = 32;             % secondary variable, used for iplot=4-6
ivar_vec = [1 2];      % set vectors to be added: uses ilev, ix, or iy as above

outfile = ('outputgraphic.eps');
outdir = ('C:\Matlab\work\images\eps_images\');

%-----------------------------------------------------------------
% Shouldn't need to change anything below here ******
%-----------------------------------------------------------------
% gnp Check if ilev1, ilev2 if in use, is inside possible levels for variable
um_level_test
if (um_error ~= 0) 
    disp('!!Terminating um_plot!!')
    return
end
% gnp Check if itime is within t:
l_time = size(t);
max_time=max(itime);
if (max_time > l_time(1))
    disp('You are asking for time outside of the t-array!')
    disp('!!Terminating um_plot!!')
    return
end
% gnp Make sure that the correct times are plotted together
% gnp T6Mean6 is plotted with the end of the mean period
um_time_test
% loop over time -------------------------------------------------------------
i = 1;
zoomey=1;  % this deals with the zooming of x_p and y_p when plotting multiple variables
for nplot = itime1
    figure(i)   
    % define colormaps, if icolormap=0 then use default (jet)
    if (icolormap == 1)
        colormap(autumn);           
        cmap1 = colormap;
        cmap2 = [cmap1(64:-1:1,1) cmap1(64:-1:1,2) cmap1(64:-1:1,3)];
        colormap(cmap2);
        brighten(0.5)
    elseif (icolormap ==2)
        colormap(gray);
        cmap1 = colormap;
        cmap2 = [cmap1(64:-1:1,1) cmap1(64:-1:1,2) cmap1(64:-1:1,3)];
        colormap(cmap2);
    end
    ilev = ilev1; 
    um_set_field;                   % field defined here  

    if zoom_plot == 1
        lsm_temp = (lsm(:,:,zooma:zoomb,zoomc:zoomd)-1).^2;
    else
        lsm_temp = (lsm-1).^2;
    end
    [aa bb] = size(lsm_temp);
    for ww=1:aa
        for qq=1:bb
            if lsm_temp(ww,qq)==0
                lsm_temp(ww,qq)=NaN;
            end
        end
    end
    if plot_over_orog == 0      % removes variable over orography if set to 0        
        plotfield = plotfield.* lsm_temp;  % removes plofield over land
    end
    
    fieldnamea = fieldname; fieldnamea2 = fieldname2; timea = time;
    fieldnamet = fieldname; fieldnamet2 = fieldname2;
    if (nplot ~= 0)                 % only plot if field found at timestep
    if (iplot == 1)                 % iplot loop generates plot 
        [C,h] = contour(dim2,dim1,plotfield,conlevs,'k');
        clabel(C,h,conlevs_labels)
        icolorbar=0;
    elseif (iplot == 2)
        [C,h] = contourf(dim2,dim1,plotfield,conlevs,'k');
        clabel(C,h,conlevs_labels)
        cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
        if (icolorbar ==1); cbh = colorbar('vert'); set(cbh,'Fontsize', 11); end;
    elseif (iplot == 3)
        pcolor(dim2,dim1,plotfield);
        cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
        shading flat;
        if (icolorbar ==1); cbh = colorbar; set(cbh,'Fontsize', 11); end;
    elseif (iplot == 4)
        [C,h] = contour(dim2,dim1,plotfield,conlevs,'k');
        clabel(C,h,conlevs_labels)
        icolorbar=0;
        hold on;
        if (ivar2 ~=0) 
            ivar_temp = ivar; ivar = ivar2; 
            nplot = itime2(i);
            ilev = ilev2;
        if (nplot ~= 0)
            um_set_field;                   % field defined here  
            [C1,h1] = contour(dim2,dim1,plotfield,conlevs,'k--');
            clabel(C1,h1,conlevs_labels)
        end
            ivar = ivar_temp;
            itime = itime1;
        end
    elseif (iplot == 5)
        [C,h] = contourf(dim2,dim1,plotfield,conlevs,'k');
%         set(h,'LineStyle','none');
%         clabel(C,h,conlevs_labels)
        cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
        if (icolorbar ==1); cbh = colorbar('vert'); set(cbh,'Fontsize', 11); end;
        hold on;
        if (ivar2 ~=0) 
            ivar_temp = ivar; ivar = ivar2; 
            nplot = itime2(i);
            ilev = ilev2;
        if (nplot ~= 0)
            um_set_field;
            [C1,h1] = contour(dim2,dim1,plotfield,conlevs,'k--');
            clabel(C1,h1,conlevs_labels)
        end
            ivar=ivar_temp;
            itime = itime1;
        end
    elseif (iplot == 6)
        pcolor(dim2,dim1,plotfield);
        cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
        shading flat;
        if (icolorbar ==1); cbh = colorbar; set(cbh,'Fontsize', 11); end;
        hold on;
        if (ivar2 ~= 0) 
            ivar_temp = ivar; ivar = ivar2; 
            nplot = itime2(i);
            ilev = ilev2;
        if (nplot ~= 0)             % only plot if field found at timestep
            um_set_field;
            [C1,h1] = contour(dim2,dim1,plotfield,conlevs,'k');
            clabel(C1,h1,conlevs_labels)
        end
            ivar=ivar_temp;
            itime = itime1;
            disp('Time for var2'); time
        end
    end
    elseif (ivar2 ~= 0 && iplot >= 4)   % ivar1 not avail. at analysis time but ivar2 is
        disp('Var1 not found at analysis time but Var2 is!')
        ivar_temp = ivar; ivar = ivar2; 
        nplot = itime2(i);
        ilev = ilev2;
        if (nplot ~= 0)             % only plot if field found at timestep
            um_set_field;
            if (iplot == 4)
                [C,h] = contour(dim2,dim1,plotfield,conlevs,'k');
                clabel(C,h,conlevs_labels)
                icolorbar=0;
            elseif (iplot == 5)        
                [C1,h1] = contour(dim2,dim1,plotfield,conlevs,'k--');
                clabel(C1,h1,conlevs_labels)
            elseif (iplot == 6)  
                [C1,h1] = contour(dim2,dim1,plotfield,conlevs,'k');
                clabel(C1,h1,conlevs_labels)
            end
        end
            itime = itime1;
            ivar=ivar_temp;
            nplot = 0;
    end

    hold on;
    set(gca,'Fontsize', 12)
    % define fieldname for title and filename
    
    if (iplot >= 4)
        fieldnameb = fieldname; fieldnameb2 = fieldname2; timeb = time;
        fieldnamet = [fieldnamea '_' fieldnameb];
        fieldnamet2 = [fieldnamea2 ' & ' fieldnameb2];
        if (timeb ~= timea) disp('Times of fields a and b are different!'); disp(timea); disp(timeb); end;
        time = t(itime)*24;
    end; 
    % add vectors
    if (ivec ~= 0) 
%         if (no_vec ~= 1)
		if (1==1)
            nplot = itime(i);
%             nplot = itime3(i);
            ilev_temp = ilev;
            ilev = ilev_vec;
            ivar_temp = ivar;
            ivar = ivar_vec(1); um_set_field;
            if plot_vect_over_orog == 0      % removes variable over orography if set to 0        
                plotfield = plotfield.* lsm_temp;  % removes plofield over land
            end       
            vec1 = plotfield;
            timec = time;        
            if (timec ~= timea); disp('Times of fields a and c (the vectors) are different!'); disp(timea); disp(timec); end;
            ivar = ivar_vec(2); um_set_field;
            if plot_vect_over_orog == 0      % removes variable over orography if set to 0        
                plotfield = plotfield.* lsm_temp;  % removes plofield over land
            end
            vec2 = plotfield;
            if (scntype ~= 1); vec2 = vec2*100; end;    % scale w velocity 
            um_add_vectors;
            ivar = ivar_temp;
            ilev = ilev_temp;
        else
            no_vec = 0;   % For all times except analyses time
        end
    end;
    % plot annotations --------------------------------------------------------
    % For x-y plots, add coastline, latitude and longitude lines, axis-labels, etc
    % Change plot aspect ratio to be equal distance in x and y directions
    if scntype==1  
        axis([min(dim2(1,:)) max(dim2(1,:)) min(dim1(:,1)) max(dim1(:,1))])
        set(gca,'PlotBoxAspectRatio', [size(dim1,2) size(dim1,1) 1]);    
        latarr = -90:5:90;
        lonarr = -180:5:180;
        set(gca,'Xtick', lonarr)
        set(gca,'Ytick', latarr)
        xlabel('longitude')
        ylabel('latitude')
        if zoom_plot==1
            [C5, h5] = contour(dim2,dim1,lsm(1,1,idim1+zooma-1,idim2+zoomc-1),[1],'k');
        else
            [C5, h5] = contour(dim2,dim1,lsm(1,1,idim1,idim2),[1],'k');     % coastline via lsm
        end
        set(h5,'Linewidth',2);
        contour(dim2,dim1,x_p(idim2,idim1)',latarr,'k');               % latitude lines    
        contour(dim2,dim1,y_p(idim2,idim1)',lonarr,'k');               % longitude lines


        if (vtype==1 || vtype==2)
            title_text = [jobname ':  ' fieldnamet2 ' at level=' num2str(ilev) ' (' z_name2 ')  ' num2str(time(i)) ' h'];
            print_text = [jobname '_' fieldnamet '_level_' num2str(ilev) '_' num2str(time(i)) 'h']
        elseif (vtype==3 | vtype==4)
            title_text = [jobname ':  ' fieldnamet2 ' at ' z_name2 '  ' num2str(time(i)) ' h'];       
            print_text = [jobname '_' fieldnamet '_' z_name '_' num2str(time(i)) 'h']
        elseif (vtype==5)
            title_text = [jobname ':  ' fieldnamet2 ' at ' num2str(time(i)) ' h'];       
            print_text = [jobname '_' fieldnamet '_' num2str(time(i)) 'h']
        end
    elseif scntype==2               % y-z plots
        axis([min(dim2(1,:)) max(dim2(1,:)) 0 height_limit])
        set(gca,'PlotBoxAspectRatio', [2 1 1]);    
        latarr = -90:5:90; 
        set(gca,'Xtick', latarr)
        set(gca,'Ytick', [1:1:height_limit]) 
        xlabel('latitude')
        if (iheight == 0)
            ylabel('hybrid height (km)')
        elseif (iheight ==1) 
            area(dim2(1,:), orography(1,1,idim2,ix)./z_scale,'Facecolor','k');      % add topography
            ylabel('height (km)')
        end; 
        if (nplot == 0); icolorbar = 0; end;
        if (icolorbar == 1); set(cbh,'position',[0.85 0.3 0.03 0.4]); end;
        title_text = [jobname ':  ' fieldnamet2 ' at ix=' num2str(ix) ' (lon=' num2str(round(mean_lon)) '^o)  ' num2str(time(i)) ' h']; 
        print_text = [jobname '_' fieldnamet '_ix_' num2str(ix) '_' num2str(iheight) '_' num2str(time(i)) 'h']
    elseif scntype==3               % x-z plots
        axis([min(dim2(1,:)) max(dim2(1,:)) 0 height_limit])
        set(gca,'PlotBoxAspectRatio', [2 1 1]);    
        lonarr = -180:5:180;
        set(gca,'Xtick', lonarr)
        set(gca,'Ytick', [1:1:height_limit]) 
        xlabel('longitude')
        if (iheight == 0)
            ylabel('hybrid height (km)')
        elseif (iheight ==1) 
            area(dim2(1,:), orography(1,1,iy,idim2)./z_scale,'Facecolor','k');      % add topography
            ylabel('height (km)')
        end; 
        if (nplot == 0); icolorbar = 0; end;
        if (icolorbar ==1); set(cbh,'position',[0.85 0.3 0.03 0.4]); end;
        title_text = [jobname ':  ' fieldnamet2 ' at iy=' num2str(iy) ' (lat=' num2str(round(mean_lat)) '^o)  ' num2str(time(i)) ' h'];       
        print_text = [jobname '_' fieldnamet '_iy_' num2str(iy) '_' num2str(iheight) '_' num2str(time(i)) 'h']
   elseif scntype==4
        axis([min(dim2(1,:)) slice_limit 0 height_limit*1000])
        set(gca,'PlotBoxAspectRatio', [2 1 1]);    
        set(gca,'Ytick', [0.5:0.5:height_limit]*1000) 
        xlabel('distance (km)')
        if (iheight == 0)
            ylabel('hybrid height (m)')
        elseif (iheight ==1)
            area(dim2(1,:), orog./z_scale,'Facecolor','k');      % add topography
            ylabel('height (m)')
        end;         
        if (nplot == 0); icolorbar = 0; end;
        if (icolorbar ==1) set(cbh,'position',[0.85 0.3 0.03 0.4]); end;
        title_text = ['noiceland_hour_' num2str(time(i))];       
        print_text = [jobname '_' fieldnamet '_interp_slice_' num2str(time(i)) 'h']
    end 
    
    % add title, print to file
    title(title_text, 'Fontsize', 12);
    if (iprint == 1) 
        print(i, '-depsc2', [outdir outfile]); 
    elseif (iprint == 2)
        print(i, '-djpeg', [outdir outfile]);
    elseif (iprint == 3)
        print(i, '-depsc2', [outdir outfile]);
    end
i=i+1;
if zoom_plot==1
    x_p=z_x_p; y_p=z_y_p;   % recover x_p and y_p
end
end     % end of nplot loop