
function [] = compare_FP_FPP_HISTORY( startyear, startmonth, endyear, endmonth, plot_time, plot_fields )

% compare monthly or monthly diurnal FP (or FPP) nc4 output from two experiments

% reichle, 23 Jan 2020

png_outpath = '/gpfsm/dnb31/rreichle/tmp_FP_FPP/plots/';

exp_path1   = '/gpfsm/dnb31/rreichle/tmp_FP_FPP/data/';
exp_path2   = '/gpfsm/dnb31/rreichle/tmp_FP_FPP/data/';

exp1        = 'f522_fp';
exp2        = 'f525land_fpp';

%file_tavg   = 'monthly';
file_tavg   = 'diurnal';

latlim = [-60, 85];

% -------------------------

if 0

  startyear   = 2019;
  startmonth  =   10;

  endyear     = 2019;
  endmonth    =   12;

  %plot_time = '530am';
  plot_time = '530pm';
  %plot_time = 'all';           % 24-hour average

  plot_fields = {'TSURF',        'tavg1_2d_lnd_Nx';    ...
		 'TSOIL1',       'tavg1_2d_lnd_Nx';    ...
		 'TSOIL2',       'tavg1_2d_lnd_Nx';    ...
		 'TSOIL3',       'tavg1_2d_lnd_Nx';    ...
		 'TSOIL4',       'tavg1_2d_lnd_Nx';    ...
		 'GWETTOP',      'tavg1_2d_lnd_Nx';    ...
		 'SFMC',         'tavg1_2d_lnd_Nx';    ...
		 'PRECTOT',      'tavg1_2d_lnd_Nx';    ...
		 'TS',           'tavg1_2d_slv_Nx';    ...
		 'T2M',          'tavg1_2d_slv_Nx';    ...
		 'T10M'          'tavg1_2d_slv_Nx';    ...
		 'QV2M',         'tavg1_2d_slv_Nx';    ...
		 'U10M',         'tavg1_2d_slv_Nx';    ...
		 'V10M',         'tavg1_2d_slv_Nx';    ...
		 'PS',           'tavg1_2d_slv_Nx';    ...
		 'FRLAND',       'const_2d_asm_Nx';    ...
		 'FRLAKE',       'const_2d_asm_Nx';    ...
		 'FRLANDICE',    'const_2d_asm_Nx';    ...
		 'FROCEAN',      'const_2d_asm_Nx';    ...
  };


  %plot_fields = {'TSURF',        'tavg1_2d_lnd_Nx';    ...
  %	       'PRECTOT',      'tavg1_2d_lnd_Nx';    ...
  %	       'GWETTOP',      'tavg1_2d_lnd_Nx';    ...
  %             };

  plot_fields = { 'PS',           'tavg1_2d_slv_Nx';    ...
		  'FRLAND',       'const_2d_asm_Nx'; };


end


% ---------------------------------------------------------------------

% load custom colormaps

mcm;


% some strings, create output directory

exp1_label = weed_out_underscores({exp1},'\_');
exp2_label = weed_out_underscores({exp2},'\_');

%out_dir =  [png_outpath, '/', exp1, '__', exp2, '/'];
out_dir = png_outpath;

eval(['!mkdir -p ', out_dir])

YYYY = num2str(startyear,  '%4.4d');
MM   = num2str(startmonth, '%2.2d');

datestringstart = [ YYYY, MM ];    


% get coordinates etc

% assume that lat/lon coordinates are the same across collections for a given experiment

% get land fraction and weights for global averages

fn1 = [exp_path1, '/', exp1, '.const_2d_asm_Nx.nc4'];
fn2 = [exp_path2, '/', exp2, '.const_2d_asm_Nx.nc4'];

lat1 = ncread(fn1,'lat');
lon1 = ncread(fn1,'lon');
lat2 = ncread(fn2,'lat');
lon2 = ncread(fn2,'lon');

if any(find(lat1-lat2)) | any(find(lon1-lon2))

  error('lats or lons in nc4 files differ.  STOPPING...')

end

FRLAND1 = ncread(fn1,'FRLAND');
FRLAND2 = ncread(fn2,'FRLAND');

% do NOT use AREA from const_2d file -- this is area of cube-sphere grid cells interpolated to lat/lon grid
%
%AREA1   = ncread(fn1,'AREA');   
%AREA2   = ncread(fn2,'AREA');

cos_lat = ones(size(lon1))*cos(lat1/180*pi)';

% weights for global data

weightsG = cos_lat/sum(cos_lat(:));

% weights for land-only data

weights1 = FRLAND1.*cos_lat;       weights1 = weights1/sum(weights1(:));
weights2 = FRLAND2.*cos_lat;       weights2 = weights2/sum(weights2(:));

weightsL = (weights1 + weights2);  weightsL = weightsL/sum(weightsL(:)); 



% ------------------------------------------------------------------------

% loop through fields and plot each 

for kk=1:size(plot_fields,1)         

  tmpvarname    = plot_fields{kk,1};
  tmpcollection = plot_fields{kk,2};

  nn = 0;

  data1    = zeros(length(lon1),length(lat1));
  data2    = zeros(length(lon2),length(lat2));

  tmpdata1 = zeros(length(lon1),length(lat1));
  tmpdata2 = zeros(length(lon2),length(lat2));

  for year=startyear:endyear

    if year==startyear

      firstmonth = startmonth;

    else

      firstmonth = 1;

    end

    if year==endyear
    
      lastmonth = endmonth;

    else

      lastmonth = 12;

    end

    for month=firstmonth:lastmonth

      YYYY = num2str(year,  '%4.4d');
      MM   = num2str(month, '%2.2d');

      datestring = [ YYYY, MM ];    

      fn1 = [exp_path1, '/', exp1, '.', tmpcollection, '.', file_tavg, '.', datestring, '.nc4'];
      fn2 = [exp_path2, '/', exp2, '.', tmpcollection, '.', file_tavg, '.', datestring, '.nc4'];

      if any(strcmp( tmpvarname, {'FRLAND','FRLAKE','FRLANDICE','FROCEAN'}))

        fn1 = [exp_path1, '/', exp1, '.', tmpcollection, '.nc4'];
        fn2 = [exp_path2, '/', exp2, '.', tmpcollection, '.nc4'];

      end


      disp(['reading ', tmpvarname, ' from ', fn1])
      disp(['reading ', tmpvarname, ' from ', fn2])
    
      tmptmpdata1 = ncread(fn1,tmpvarname);
      tmptmpdata2 = ncread(fn2,tmpvarname);

      % remove GEOS no-data-values

      tmptmpdata1( tmptmpdata1 > .9999e15  ) = NaN;
      tmptmpdata2( tmptmpdata2 > .9999e15  ) = NaN;

      % change units of precip
   
      if strcmp(tmpvarname,'PRECTOT')

        tmptmpdata1 = tmptmpdata1*86400;
        tmptmpdata2 = tmptmpdata2*86400;

      end

      if any(strcmp( tmpvarname, {'FRLAND','FRLAKE','FRLANDICE','FROCEAN'}))

	  tmpdata1 = tmptmpdata1;
          tmpdata2 = tmptmpdata2;

      else

        % compute daily average or extract data for a given local solar time
  
        % "diurnal" files contain 24 times: 0:30z, 1:30z, 2:30z, ..., 23:30z
        %                            index: 1      2      3           24

        switch plot_time
  
          case '430am'
  
            offset = 4;
  
          case '530am'
  
            offset = 5;

          case '630am'
  
            offset = 6;
  
          case '730am'
  
            offset = 7;
  

          case '430pm'
  
            offset = -8;

          case '530pm'
  
            offset = -7;

          case '630pm'
  
            offset = -6;

          case '730pm'
  
            offset = -5;
  
        end
  
  
        switch plot_time
  
          case {'430am','530am','630am','730am',...
                '430pm','530pm','630pm','730pm'}
  
            tmpind = mod( round((360-lon1)/15)+offset, 24) * ones(size(lat1') ) + 1; 
  
            for ii=1:length(lon1)
  	      for jj=1:length(lat1)
  	      
  	        tmpdata1(ii,jj) = tmptmpdata1( ii, jj, tmpind(ii,jj) );
  	        tmpdata2(ii,jj) = tmptmpdata2( ii, jj, tmpind(ii,jj) );
  
              end
  	    end
  
          case 'all'
  
	    tmpdata1 = mean(tmptmpdata1,3);
            tmpdata2 = mean(tmptmpdata2,3);

          otherwise
  
            error('unknown "plot_time", STOPPING...')
  
        end

      end

      % sum up for average

      data1 = data1 + tmpdata1;
      data2 = data2 + tmpdata2;
  
      nn = nn + 1; 
   
    end

  end 

  % finalize average

  data1 = data1 / nn; 
  data2 = data2 / nn;

  % compute difference field

  datadiff = data2 - data1;   

  % assemble datestring for averaging period

  datestring = [datestringstart, '-', datestring];

  % read units

  tmpunits    = ncreadatt(fn1, tmpvarname,'units');
  tmplongname = ncreadatt(fn1, tmpvarname,'long_name');

  if strcmp(tmpvarname,'PRECTOT'), tmpunits = 'mm d-1'; end

  % get stats

  datamax  = max( [data1(:); data2(:)] );
  datamin  = min( [data1(:); data2(:)] );
  
%  if findstr(tmpcollection,'lnd')
%
%    aavediff  = nansum1d(  datadiff(:).*weightsL(:) );
%
%  else
%
%    aavediff  = nansum1d(  datadiff(:).*weightsG(:) );
%
%  end
%
%  rmsd        = nanstd1d(  datadiff(:) );


  if contains(tmpcollection,'lnd')

    [aavediff, rmsd]  = wmeannan(  datadiff(:), weightsL(:) );

  else

    [aavediff, rmsd]  = wmeannan(  datadiff(:), weightsG(:) );

  end


  maxabsdiff  = max( abs(  datadiff(:)));
    

  tmplongname = weed_out_underscores({tmplongname}, ' ');

  tmp_title   = [ plot_time, ': ',  tmpvarname, ' (', tmplongname{1}, ')  [', tmpunits, '],  ', datestring ];

  tmp_titlediff = {['Top minus Middle']; 
		   ['aave =',    num2str(aavediff),   ',   ',                      ... 
	            'RMSD=',     num2str(rmsd),       ',   ',                      ... 
                    'MaxAbs=',   num2str(maxabsdiff)                               ...
                   ]};
  

  % colormap for fields

  tmpcm    = parula(64);
 
  tmpcminv = tmpcm(end:-1:1,:);

  % colormap for diff maps

  %cmax = maxabsdiff;
  cmax = 3*rmsd;

  if cmax==0,  cmax = 1e-10;  end
    
  tmpcmdiff = mcm_BuDRd_13;
  
  tmpcmdiff(7,:) = [.5 .5 .5];
   
  tmpcmdiffinv = tmpcmdiff(end:-1:1,:);

  % pick appropriate colormap

  if any(strcmp( tmpvarname, {'GWETTOP','SFMC','PRECTOT','QV2M'}))

    thiscm     = tmpcminv;
    thiscmdiff = tmpcmdiffinv;

  else 

    thiscm     = tmpcm;
    thiscmdiff = tmpcmdiff;

  end


  % ---------------------------------------------

  % scatter plot

  plot_scatter = 0;

  if plot_scatter

    hf = figure;
  
    plot( data2(:), data1(:), '.')
  
    title(tmp_title);

    xlabel( exp2_label )
    ylabel( exp1_label )
    
    % print to file
  
    drawnow;
  
    fname_out = [out_dir, '/scatter_', datestring, '_', tmpvarname, '_', plot_time, '.png'];
  
    print(fname_out,'-dpng')
  
    close(hf)

  end
  
  % ---------------------------------------------

  % map plots

  plot_maps = 1;

  if plot_maps

    hf = figure;  orient tall;  colormap(thiscm);

    subplot(3,1,1)

    pcolor( lon2, lat2, data2' )
  
    shading('flat'), colorbar

    if findstr(tmpcollection,'lnd'), set(gca,'YLim',latlim); end

    title( {tmp_title;  exp2_label{1} } )

    caxis([datamin, datamax])  



    subplot(3,1,2)

    pcolor( lon1, lat1, data1' )      

    shading('flat'), colorbar

    if findstr(tmpcollection,'lnd'), set(gca,'YLim',latlim); end

    title( exp1_label{1} )

    caxis([datamin, datamax])  



    subplot(3,1,3)

    pcolor( lon1, lat1, datadiff' )
  
    shading('flat'), colorbar

    if findstr(tmpcollection,'lnd'), set(gca,'YLim',latlim); end
  
    title( tmp_titlediff )
  
    caxis([-cmax cmax])  

    colormap(gca, thiscmdiff)



    drawnow;
  
    fname_out = [out_dir, '/map_', datestring, '_', tmpvarname, '_', plot_time, '.png'];
  
    disp(['printing ', fname_out])

    print(fname_out,'-dpng')
  
    close(hf)

  end
  
end


  
% =================== EOF ==================================


