% um_set_field.m
% sets plotfield, fieldname and dim1, dim2, plus vtype and conlevs
% called from um_plot

%-----------------------------------------------------------------
% set field, fieldnames, vertical coordinate and contour intervals 
%-----------------------------------------------------------------
% 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.......
%-----------------------------------------------------------------
clear field1 plotfield conlevs conlevs_labels pf d1 d2 XI YI ZI column temp_orog slice_z zd slice
time =t(itime)*24
if ivar==1           
    field1 = u; fieldname = 'u'; fieldname2= 'u velocity (m s^{-1})';
    vtype = 1; 
    conlevs = -50:5:50; conlevs_labels = -50:10:50;   
elseif ivar == 2
    field1 = v; fieldname = 'v'; fieldname2 = 'v velocity (m s^{-1})';
    vtype = 1; 
    conlevs = -30:5:30; conlevs_labels = -30:10:30;   
elseif ivar == 3
    field1 = theta; fieldname = 'theta'; fieldname2 = 'potential temperature (K)';
    vtype = 2; 
    conlevs = 264:2:300; conlevs_labels = 264:4:300;
elseif ivar == 4
    field1 = q*1000; fieldname = 'q'; fieldname2 = 'specific humidity (g kg^{-1})';
    vtype = 2; 
    conlevs = 0:0.2:5; conlevs_labels = 0:0.4:5;
elseif ivar == 5
    field1 = w; fieldname = 'w'; fieldname2 = 'w velocity (m s^{-1})';
    vtype = 2;
    conlevs = -0.7:0.05:0.3; conlevs_labels = -0.7:0.05:0.3; 
elseif ivar == 6    
    field1 = pv*1e6; fieldname = 'pv'; fieldname2 = 'potential vorticity (PVU)';
    vtype = 1;
    conlevs = [-6 -4 -2 -1 0 1 2 4 6 8]; conlevs_labels = conlevs;
elseif ivar == 7    
    field1 = qcl*1000; fieldname = 'qcl'; fieldname2 = 'cloud liquid water (g kg^{-1})';
    vtype = 2; 
    conlevs = -0.2:0.005:0.2; conlevs_labels = -0.2:0.01:0.2;
elseif ivar == 8
    field1 = qcf*1000; fieldname = 'qcf'; fieldname2 = 'cloud ice content (g kg^{-1})';
    vtype = 2;
    conlevs = -0.2:0.05:0.6; conlevs_labels = -0.2:0.05:0.6;
elseif ivar == 9
    field1 = pv_theta*1e6; fieldname = 'pv_theta'; fieldname2 = 'potential vorticity (PVU)';
    vtype = 3;
    conlevs = [-8 -6 -4 -2 -1 0 1 2 4 6 8]; conlevs_labels = conlevs;
elseif ivar == 10
    field1 = gph; fieldname = 'gph'; fieldname2 = 'geopotential height (m)';
    vtype = 4;
    conlevs = 5000:10:5500; conlevs_labels = 5000:20:5500;
elseif ivar == 11
    field1 = lsm; fieldname = 'lsm'; fieldname2 = 'land sea mask';
    vtype = 5; %ilev=1;
    conlevs = 0.5:0.5:1; conlevs_labels = NaN;      
elseif ivar == 12
    field1 = orography; fieldname = 'orography'; fieldname2 = 'orography (m)';
    vtype = 5; %ilev=1;
    conlevs = 0:250:4000; conlevs_labels = 0:500:4000;  
elseif ivar == 13
    field1 = temp_sfc; fieldname = 'temp_sfc'; fieldname2 = 'surface temperature (K)';
    vtype = 5; %ilev=1;
    conlevs = 220:2:300; conlevs_labels = 220:2:300;    
elseif ivar == 14
    field1 = mslp/100; fieldname = 'mslp'; fieldname2 = 'msl pressure (mb)';
    vtype = 5; %ilev=1;
    conlevs = 940:4:1040; conlevs_labels = 904:8:1040;
elseif ivar == 15
    field1 = temp_1pt5m; fieldname = 'temp_1pt5m'; fieldname2 = 'temperature at 1.5 m (K)';
    vtype = 5; %ilev=1;
    conlevs = 220:2:286; conlevs_labels = 220:2:286;   
elseif ivar == 16
    field1 = q_1pt5m*1000; fieldname = 'q_1pt5m'; fieldname2 = 'specific humidity at 1.5 m (g kg^{-1})';
    vtype = 5; %ilev=1;
    conlevs = 0:0.5:10; conlevs_labels = 0:1:10;   
elseif ivar == 17
    field1 = u_10m; fieldname = 'u_10m'; fieldname2= 'u_{10m} (m s^{-1})';
    vtype = 5; %ilev=1;
    conlevs = -50:5:50; conlevs_labels = -50:10:50;   
elseif ivar == 18
    field1 = v_10m; fieldname = 'v_10m'; fieldname2= 'v_{10m} (m s^{-1})';
    vtype = 5; %ilev=1;
    conlevs = -50:5:50; conlevs_labels = -50:10:50;   
elseif ivar == 19
    field1 = wind_speed_10m; fieldname = 'wind_speed_10m'; fieldname2= 'wind speed at 10 m (m s^{-1})';
    vtype = 5; %ilev=1;
    conlevs = 0:2:35; conlevs_labels = 0:4:35;  
elseif ivar == 20
    field1 = low_cloud; fieldname = 'low_cloud'; fieldname2= 'low cloud (fraction)';
    vtype = 5; %ilev=1;
    conlevs = 0:0.2:1; conlevs_labels = 0:0.4:1;  
elseif ivar == 21
    field1 = med_cloud; fieldname = 'med_cloud'; fieldname2= 'medium cloud (fraction)';
    vtype = 5; %ilev=1;
    conlevs = 0:0.2:1; conlevs_labels = 0:0.4:1;  
elseif ivar == 22
    field1 = high_cloud; fieldname = 'high_cloud'; fieldname2= 'high cloud (fraction)';
    vtype = 5; %ilev=1;
    conlevs = 0:0.2:1; conlevs_labels = 0:0.4:1;  
elseif ivar == 23
    field1 = shf_sfc; fieldname = 'shf_sfc'; fieldname2 = 'mean surface sensible heat flux (W m^{-2})';
    vtype = 5; %ilev=1;
    conlevs = -200:50:300; conlevs_labels = -200:50:300; 
elseif ivar == 24
    field1 = lhf_sfc; fieldname = 'lhf_sfc'; fieldname2 = 'mean surface latent heat flux (W m^{-2})';
    vtype = 5; %ilev=1;
    conlevs = -200:50:300; conlevs_labels = -200:50:300;
elseif ivar == 25
    field1 = taux_sfc; fieldname = 'taux_sfc'; fieldname2 = 'mean surface u-momentum flux (N m^{-2})';
    vtype = 5; %ilev=1;
    conlevs = 0:0.1:1; conlevs_labels = 0:0.1:1;
elseif ivar == 26
    field1 = tau_sfc; fieldname = 'tauy_sfc'; fieldname2 = 'mean surface v-momentum flux (N m^{-2})';
    vtype = 5; %ilev=1;
    conlevs = 0:0.1:1; conlevs_labels = 0:0.1:1;
elseif ivar == 27
    field1 = tau_sfc; fieldname = 'tau_sfc'; fieldname2 = 'mean surface momentum flux (N m^{-2})';
    vtype = 5; %ilev=1;
    conlevs = 0:0.2:3; conlevs_labels = 0:0.4:3;  
elseif ivar == 28
    field1 = lsrain; fieldname = 'lsrain'; fieldname2 = 'accumulated large scale rain (mm)';
    vtype = 5; %ilev=1;
    conlevs = 1:0.5:3; conlevs_labels = conlevs;  
elseif ivar == 29
    field1 = lssnow; fieldname = 'lssnow'; fieldname2 = 'accumulated large scale snow (mm w.e.)';
    vtype = 5; %ilev=1;
    conlevs = 1:1:10; conlevs_labels = conlevs;  
elseif ivar == 30
    field1 = cvrain; fieldname = 'cvrain'; fieldname2 = 'accumulated convective rain (mm)';
    vtype = 5; %ilev=1;
    conlevs = 1:1:10; conlevs_labels = conlevs;  
elseif ivar == 31
    field1 = cvsnow; fieldname = 'cvsnow'; fieldname2 = 'accumulated convective snow (mm w.e.)';
    vtype = 5; %ilev=1;
    conlevs = 1:1:10; conlevs_labels = conlevs;  
elseif ivar == 32
    field1 = wind_speed; fieldname = 'wind speed'; fieldname2 = 'wind speed (m s^{-1})';
    vtype = 1; time =t(itime)*24;
    conlevs = 0:2:50; conlevs_labels = 0:4:50; 
elseif ivar == 33
    field1 = p_levels1./100; fieldname = 'pressure'; fieldname2 = 'pressure (mb)';
    vtype = 1; 
    conlevs = 600:2:1050; conlevs_labels = 600:4:1050;   
elseif ivar == 34
    field1 = wind_dir; fieldname = 'wind direction'; fieldname2 = 'wind direction (^o)';
    vtype = 1; time =t(itime)*24;
    conlevs = 0:10:360; conlevs_labels = 0:20:360; 
elseif ivar == 35
    field1 = ice_conc.*100; fieldname = 'ice conc'; fieldname2 = 'ice concentration (%)';
    vtype = 5; %ilev=1;
    conlevs = 0:10:100; conlevs_labels = 0:20:100;    
end;

% define vertical coordinate 
if vtype==1 
    if (iheight == 0)           % dim1 = hybrid height
        zdim=hheight_rho; 
    elseif (iheight == 1)       % dim1 = geometric height
        zdim = z_rho(1,:,:,:);
    end; 
    z_name = num2str(round(hheight_rho(ilev))); z_name2 =[z_name ' m']; 
elseif vtype==2 
    if (iheight == 0)           % dim1 = hybrid height
        zdim=hheight_theta;    
    elseif (iheight == 1)       % dim1 = geometric height
        zdim = z_theta(1,:,:,:);
    end;
    z_name = num2str(round(hheight_theta(ilev))); z_name2 =[z_name ' m']; 
elseif vtype==3 
    zdim=theta_levels; z_name = num2str(round(zdim(ilev))); z_name2 =[z_name ' K']; 
elseif vtype==4 
    zdim=p_levels; z_name = num2str(round(zdim(ilev))); z_name2 =[z_name ' mb']; 
elseif vtype==5
    zdim= 0;
end

% reduces variable being plotted to the zoomed area 
% and does the same for the x_p and y_p coordinates
if zoom_plot==1 & zoomey==1
    z_x_p=x_p; z_y_p=y_p;   % store x_p and y_p
    field1 = field1(:,:,zooma:zoomb,zoomc:zoomd); 
    x_p=x_p(zoomc:zoomd,zooma:zoomb); % **************** added to allow for zoom
    y_p=y_p(zoomc:zoomd,zooma:zoomb);
    zoomey=zoomey+1;
elseif zoom_plot==1 & zoomey==2
    field1 = field1(:,:,zooma:zoomb,zoomc:zoomd); 
end

% determines size of domain from plotting variable
idimxmin = 1; idimymin =1 ;
sizedim = size(field1);
if ivar == 12
    Nx = sizedim(2); Ny = sizedim(1);    %  used for orography
else
    Nt = sizedim(1); Nz = sizedim(2); Nx = sizedim(4); Ny = sizedim(3);
end

% define domain of variable: x-y, y-z, x-z, distance-z (scntype 1,2,3,4)
% i.e. defined plotfield, dim1 and dim2
% dim1 and dim2 now defined as matrices for all scntypes
% size Ny by Nx, Nz by Ny, Nz by Nx or Nz by distance
z_scale = 1000;
if scntype == 1               % creates x-y plot
    idim2 = 1:Nx; idim1 = 1:Ny; 
    dim2 = x_p(idim2,idimymin); dim1 = y_p(idimxmin,idim1);
    dim2 = repmat(dim2,1,Ny)'; dim1 = repmat(dim1,Nx,1)';
    nplot,ilev
    if (nplot ~= 0)
        if ivar == 12
            plotfield(idim1,idim2) = field1;        % used for orography 
        else
            plotfield(idim1,idim2) = field1(nplot,ilev,idim1,idim2);
        end
    end
elseif scntype == 2            % creates y-z plot
    idim2 = 1:Ny; idim1 = 1:Nz;
    dim2 = y_p(ix,idim2);
    mean_lon = mean(x_p(ix,:)); disp(['Mean longitude = ' num2str(mean_lon)]);
    if (iheight ==0)
        dim1 = zdim./z_scale;
        dim1 = repmat(dim1,1,Ny); 
    elseif (iheight == 1)
        dim1 = zdim(idim1,idim2,ix)./z_scale;
    end;  
    dim2 = repmat(dim2,Nz,1);
    if (nplot ~= 0)
        plotfield(idim1,idim2) = field1(nplot,idim1,idim2,ix);
    end
elseif scntype == 3              % creates x-z plot
    idim2 = 1:Nx; idim1 = 1:Nz; 
    dim2 = x_p(idim2,iy); 
    mean_lat = mean(y_p(:,iy)); disp(['Mean latitude = ' num2str(mean_lat)]);
    if (iheight ==0)
        dim1 = zdim./z_scale;
        dim1 = repmat(dim1,1,Nx); 
    elseif (iheight == 1)
        dim1 = zdim(idim1,iy,idim2)./z_scale;
        dim1 = reshape(dim1,Nz,Nx);
    end;        
    dim2 = repmat(dim2,1,Nz)';    
    if (nplot ~=0)
        plotfield(idim1,idim2) = field1(nplot,idim1,iy,idim2);
    end
elseif scntype == 4           % creates and interpolated vertical cross section between two points (AX,AY) and (BX,BY)
    temp_orog(:,:) = orography(1,1,:,:);
    temp_zdim = zdim*1000;
    dlength=12;     % ******** domain grid resolution = 12km
	if AX>BX            % *********** extracts vertical column of data between the two points
        column=field1(:,:,:,BX-10:AX+10);
        z_column=temp_zdim(:,:,BX-10:AX+10);
        temp_orog=temp_orog(:,BX-10:AX+10);
        x1=11; x2=(AX-BX)+11;
	else
        column=field1(:,:,:,AX-10:BX+10);
        z_column=temp_zdim(:,:,AX-10:BX+10);
        temp_orog=temp_orog(:,AX-10:BX+10);
        x1=11; x2=(BX-AX)+11;
	end
	if AY>BY
        column=column(:,:,BY-10:AY+10,:);
        z_column=z_column(:,BY-10:AY+10,:);
        temp_orog=temp_orog(BY-10:AY+10,:);
        y1=11; y2=(AY-BY)+11;
	else
        column=column(:,:,AY-10:BY+10,:);
        z_column=z_column(:,AY-10:BY+10,:);
        temp_orog=temp_orog(AY-10:BY+10,:);
        y1=11; y2=(BY-AY)+11;
	end
    
    % If it is a NW-SE slice then rotation needed in order to get the data on
    % the diagonal
    if (AY>BY && AX<BX) || (AY<BY && AX>BX)
        rot=-1;
        disp('Slices: rotating the matrix:)')
    else
        rot=0;
    end
    
    % seting up the horizontal grid for 2-D interpolation
    % interpolates from the model domain onto a grid between the two points
    % extracts the diagonal as first level of vertical cross section
	[d2, d1] = meshgrid(x1-10:x2+10,y1-10:y2+10); 
    distance=dlength*sqrt((x2-x1)^2+(y2-y1)^2);
	ratio=((x2-x1)/(y2-y1));
    dy=sqrt((dl^2)/(ratio^2+1)); % ******************* n_x is not a whole number, but this is correct
	dx=sqrt(dl^2-dy^2);
	n_x=((x2-x1)*dlength)/dx;
	n_y=((y2-y1)*dlength)/dy;
	x2=x1+(ceil(n_x)*dx)/dlength;
	y2=y1+(ceil(n_y)*dy)/dlength;
    [XI,YI] = meshgrid(x1:dx/dlength:x2,y1:dy/dlength:y2);
  
    ZI = interp2(d2,d1,temp_orog,XI,YI);
    ZI = rot90(ZI,rot);
    orog=diag(ZI);

%     ****************** used for checking interpolation scheme ************************    
%     pf(:,:)=column(3,1,:,:);  %***************** usually commented ******************
%     figure(1)
%     pcolor(d2,d1,pf)
%     cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
%     shading flat;
%     figure(2)
%     ZI = interp2(d2,d1,pf,XI,YI);
%     pcolor(XI,YI,ZI)
%     cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
%     shading flat;
%     axis([x1 x2 y1 y2 -5 10])
%     view([0 90])
%     clear pf              %***************** usually commented ******************
%     *******************************************************************************

    Nslice = ceil(n_x)+1;
    idim2 = 1:Nslice; idim1 = 1:Nz;
    dim2 = (idim2-1)'.*dl;           % -1 so that the first point is at 0km    
    if (iheight ==0)
        dim1 = zdim./z_scale;
        dim1 = repmat(dim1,1,Nslice);
    end
    
    % repeats the above interpolation for the remaining levels
    % creating a vertical cross section in slice
    for int_loop=1:Nz 
        pf(:,:)=column(nplot,int_loop,:,:);
        zd(:,:)=z_column(int_loop,:,:);
        ZI = interp2(d2,d1,pf,XI,YI,'spline');
        ZZ = interp2(d2,d1,zd,XI,YI,'spline');
        ZI=rot90(ZI,rot);
        ZZ=rot90(ZZ,rot);
        slice(:,int_loop)=diag(ZI);
        slice_z(:,int_loop)=diag(ZZ)./z_scale;
    end
    slice=slice';
    slice_z=slice_z';

%     ********************** used for checking interpolation scheme ***************    
%     figure(3)             %***************** usually commented ******************
%     pcolor(slice)
%     cmin=min(conlevs);cmax=max(conlevs); caxis([cmin cmax]);caxis manual;
%     shading flat;        %***************** usually commented ******************
%     ****************************************************************************** 
    
    % sets up dim1 and dim2 for plotting the cross section
    % and sets up plotfield as the cross section to be plotted
    dim1 = slice_z;
    dim2 = repmat(dim2,1,Nz)';    
	if ivar==32
        [aa bb]=size(slice);
        for temp_i=1:aa
            for temp_j=1:bb
                if slice(temp_i,temp_j)<0
                    slice(temp_i,temp_j)=0;
                end
            end
        end
    end
    plotfield = slice;
    
end