% Compare models with flux analyses and the NOAA PSD observations
% in the eastern Pacific along 20S in October. This program compares two
% models (one in each column of plots) with the observations, but could be
% extended to more models.
%
% Uses mexnc and the old netcdf toolbox.
% 
% The following data sets are required:
% WHOI OAFlux gridded flux climatology (Yu and Weller 2007)
%   .\data\lhsh_oaflux_clim_1984-2002.nc (turbulent)
%   .\data\lwsw_isccp_clim_1984-2002.nc (radiative from ISCCP FD)
% NCAR CORE gridded flux climatology (Large and Yeager 2008).
%   .\data\core_clim_1984-2004.nc
% PSD ship observations averaged for October along 20S.
%   .\data\psdobs.mat
% CMIP3 model output. It is arranged in these directories:
%   [datadir]\flux\[prefix][model][suffix][var]*.nc
%            \radiation\[model]\[var]_*.nc
%            \sftlf\[model]_*.nc
%
% (c) Simon de Szoeke
%     2008 April 2
%     revised for distribution 2009 Jan 6

latlim=[-25 -15];
lonlim=[-110 -70]+360;
imonth=10; % October

datadir='C:\Data\ipcc\20c3m\';  % for model files
prefix='pcmdi.ipcc4.';
suffix='.20c3m.run*.monthly.';

numncobs=2; % number of netcdf non-model files to be read
numncmodels=2;
model={'whoi'
       'core'
       'gfdl_cm2_0'
       'gfdl_cm2_1'};
name=lower({'WHOI/ISCCP'
            'CORE'
            'GFDL CM2.0'
            'GFDL CM2.1'});
var={'hfls'
     'hfss'
     'rlds'
     'rlus'
     'rsds'
     'rsus' 
     'sftlf'};

clf
% WHOI, Yu and Weller 2007
nmod=1; % used as a subscript to the netcdf objects
nc{nmod,1}=netcdf('.\data\lhsh_oaflux_clim_1984-2002.nc');
lat=nc{nmod,1}{'lat'}(:);
lon=nc{nmod,1}{'lon'}(:);
filat=find(lat>=latlim(1) & lat<=latlim(2));
filon=find(lon>=lonlim(1) & lon<=lonlim(2));
lh=nc{nmod,1}{'lhtfl'}(:,filat,filon);
sh=nc{nmod,1}{'shtfl'}(:,filat,filon);
lh(lh>1e4)=NaN;
sh(sh>1e4)=NaN;
nc{nmod,2}=netcdf('.\data\lwsw_isccp_clim_1984-2002.nc');
lw=nc{nmod,2}{'nlwrs'}(:,filat,filon);
sw=nc{nmod,2}{'nswrs'}(:,filat,filon);
lw(lw>1e4)=NaN;
sw(sw>1e4)=NaN;
i20s=find(abs(lat(filat)--20)<=0.5+eps);
lh20s=squeeze(mean(lh(:,i20s,:),2));
sh20s=squeeze(mean(sh(:,i20s,:),2));
lw20s=squeeze(mean(lw(:,i20s,:),2));
sw20s=squeeze(mean(sw(:,i20s,:),2));
net20s=lh20s+sh20s+lw20s+sw20s;
% plot WHOI in panels for all models
for ncol=1:numncmodels
    sax(ncol,1)=subplot(5,2,0+ncol,'align');
    hline(ncol,1)=plot(360-lon(filon),lh20s(imonth,:),'r-'); % Oct.
    hold on
    sax(ncol,2)=subplot(5,2,2+ncol,'align');
    hline(ncol,2)=plot(360-lon(filon),sh20s(imonth,:),'r-'); % Oct.
    hold on
    sax(ncol,3)=subplot(5,2,4+ncol,'align');
    hline(ncol,3)=plot(360-lon(filon),lw20s(imonth,:),'r-'); % Oct.
    hold on
    sax(ncol,4)=subplot(5,2,6+ncol,'align');
    hline(ncol,4)=plot(360-lon(filon),sw20s(imonth,:),'r-'); % Oct.
    hold on
    sax(ncol,5)=subplot(5,2,8+ncol,'align');
    hline(ncol,5)=plot(360-lon(filon),net20s(imonth,:),'r-'); % Oct.
    hold on
end % plot WHOI

% CORE, Large and Yeager 2008
nmod=2;
nc{nmod,1}=netcdf('.\data\core_clim_1984-2004.nc');
lat=nc{nmod,1}{'lat'}(:);
lon=nc{nmod,1}{'lon'}(:);
filat=find(lat>=latlim(1) & lat<=latlim(2));
filon=find(lon>=lonlim(1) & lon<=lonlim(2));
LWDN_a2c=nc{nmod,1}{'Q_lwdn'}(:,filat,filon);
SWNET_a2c=nc{nmod,1}{'Q_swnet'}(:,filat,filon);
LWUP_c2a=nc{nmod,1}{'Q_lwup'}(:,filat,filon);
lat_c2a=nc{nmod,1}{'Q_lat'}(:,filat,filon);
sens_c2a=nc{nmod,1}{'Q_sen'}(:,filat,filon);
i20s=find(abs(lat(filat)--20)<=0.5+eps);
lh20s=squeeze(mean(lat_c2a(:,i20s,:),2));
sh20s=squeeze(mean(sens_c2a(:,i20s,:),2));
lw20s=squeeze(mean(LWDN_a2c(:,i20s,:)+LWUP_c2a(:,i20s,:),2));
sw20s=squeeze(mean(SWNET_a2c(:,i20s,:),2));
net20s=lh20s+sh20s+lw20s+sw20s;
% plot CORE to all panels
for ncol=1:numncmodels;
    axes(sax(ncol,1));
    hline(ncol,1)=plot(360-lon(filon),lh20s(imonth,:),'k-'); % Oct.
    hold on
    axes(sax(ncol,2))
    hline(ncol,2)=plot(360-lon(filon),sh20s(imonth,:),'k-'); % Oct.
    hold on
    axes(sax(ncol,3));
    hline(ncol,3)=plot(360-lon(filon),lw20s(imonth,:),'k-'); % Oct.
    hold on
    axes(sax(ncol,4));
    hline(ncol,4)=plot(360-lon(filon),sw20s(imonth,:),'k-'); % Oct.
    hold on
    axes(sax(ncol,5));
    hline(ncol,5)=plot(360-lon(filon),net20s(imonth,:),'k-'); % Oct.
    hold on
end % plot CORE

for nmod=(1:numncmodels)+numncobs; % load GFDL CM2.0 and 2.1 CMIP3 models
    for nvar=1:2
        file=dir([datadir 'flux\' prefix model{nmod} suffix var{nvar} '*.nc']);
        nc{nmod,nvar}=netcdf([datadir 'flux\' file.name]);
    end
    for nvar=3:6
        file=dir([datadir 'radiation\' model{nmod} '\' var{nvar} '_*.nc']);
        nc{nmod,nvar}=netcdf([datadir 'radiation\' model{nmod} '\' file.name]);
    end
    nvar=7;
        file=dir([datadir 'sftlf\' model{nmod} '_*.nc']);
        nc{nmod,nvar}=netcdf([datadir 'sftlf\' file.name]);
end

% CMIP plot: GFDL CM2.0 and 2.1
for nmod=(1:numncmodels)+numncobs;
    lat=nc{nmod,1}{'lat'}(:);
    lon=nc{nmod,1}{'lon'}(:);
    filat=find(lat>=latlim(1) & lat<=latlim(2));
    filon=find(lon>=lonlim(1) & lon<=lonlim(2));
    lh=-nc{nmod,1}{'hfls'}(:,filat,filon);
    sh=-nc{nmod,2}{'hfss'}(:,filat,filon);
    lwd=nc{nmod,3}{'rlds'}(:,filat,filon);
    lwu=nc{nmod,4}{'rlus'}(:,filat,filon);
    swd=nc{nmod,5}{'rsds'}(:,filat,filon);
    swu=nc{nmod,6}{'rsus'}(:,filat,filon);
    sftlf=nc{nmod,7}{'sftlf'}(filat,filon);
    lh(repmat(shiftdim(sftlf,-1),[size(lh,1) 1 1]) | abs(lh)>1e4)=NaN;
    sh(repmat(shiftdim(sftlf,-1),[size(sh,1) 1 1]) | abs(sh)>1e4)=NaN;
    lw=lwd-lwu;
    sw=swd-swu;
    lw(repmat(shiftdim(sftlf,-1),[size(lw,1) 1 1]) | abs(lwd)>1e4 | abs(lwu)>1e4)=NaN;
    sw(repmat(shiftdim(sftlf,-1),[size(sw,1) 1 1]) | abs(swd)>1e4 | abs(swu)>1e4)=NaN;
    
    i20s=find(abs(lat(filat)--20)<=1.0+eps);
    if isempty(i20s)
        [junk,i20s]=min(abs(lat(filat)--20));
    end
    lh20s=NaN+zeros(12,length(filon));
    sh20s=lh20s; lw20s=lh20s; sw20s=lh20s;
    for month=1:12 % average climatology
        lh20s(month,:)=squeeze(nanmean(nanmean(lh(month:12:end,i20s,:)),2));
        sh20s(month,:)=squeeze(nanmean(nanmean(sh(month:12:end,i20s,:)),2));
        lw20s(month,:)=squeeze(nanmean(nanmean(lw(month:12:end,i20s,:)),2));
        sw20s(month,:)=squeeze(nanmean(nanmean(sw(month:12:end,i20s,:)),2));
    end
    net20s=lh20s+sh20s+lw20s+sw20s;

    ncol=nmod-numncobs; 
    axes(sax(ncol,1));
    hline(nmod,1)=plot(360-lon(filon),lh20s(imonth,:),'k-','linewidth',1.5);
    hold on
    axes(sax(ncol,2))
    hline(nmod,2)=plot(360-lon(filon),sh20s(imonth,:),'k-','linewidth',1.5);
    hold on
    axes(sax(ncol,3))
    hline(nmod,3)=plot(360-lon(filon),lw20s(imonth,:),'k-','linewidth',1.5);
    hold on
    axes(sax(ncol,4));
    hline(nmod,4)=plot(360-lon(filon),sw20s(imonth,:),'k-','linewidth',1.5);
    hold on
    axes(sax(ncol,5));
    hline(nmod,5)=plot(360-lon(filon),net20s(imonth,:),'k-','linewidth',1.5);
    hold on
end % END CMIP

% PSD ship observations
psdobs=load('.\data\psdobs.mat');
sqtn=sqrt(length(psdobs.yroct));
for ncol=1:numncmodels
    errorbar(sax(ncol,1),-psdobs.lonctr,psdobs.lhmean,psdobs.lhstd/sqtn,'b.-')
    errorbar(sax(ncol,2),-psdobs.lonctr,psdobs.shmean,psdobs.shstd/sqtn,'b.-')
    errorbar(sax(ncol,3),-psdobs.lonctr,psdobs.lwmean,psdobs.lwstd/sqtn,'b.-')
    errorbar(sax(ncol,4),-psdobs.lonctr,psdobs.swmean,psdobs.swstd/sqtn,'b.-')
    errorbar(sax(ncol,5),-psdobs.lonctr,psdobs.net,psdobs.swstd/sqtn,'b.-')
end

%WHOI Buoy obs, read off Al Plueddemann's plot
wlh=-100;
wsh=-5;
wlw=-35;
wsw= 210;
wnet=60;

for ncol=1:numncmodels
    plot(sax(ncol,1),85.5,wlh,'mo')
    plot(sax(ncol,2),85.5,wsh,'mo')
    plot(sax(ncol,3),85.5,wlw,'mo')
    plot(sax(ncol,4),85.5,wsw,'mo')
    plot(sax(ncol,5),85.5,wnet,'mo')
end

% finish plot
for ncol=1:numncmodels % column titles
    nmod=ncol+numncobs;
    set(get(sax(ncol,1),'title'),'string',name{nmod});
end
set(get(sax(1,1),'ylabel'),'string','latent');
set(get(sax(1,2),'ylabel'),'string','sensible');
set(get(sax(1,3),'ylabel'),'string','longwave');
set(get(sax(1,4),'ylabel'),'string','solar');
set(get(sax(1,5),'ylabel'),'string','net');
set(get(sax(1,5),'xlabel'),'string','west longitude')
%set(hline(:),'clipping','off')
set(sax(:),'xdir','reverse','xlim',[70 110],'xtick',70:15:110,'tickdir','out')
set(sax(:,1),'ylim',[-140 -20],'ytick',-160:40:-20);
set(sax(:,2),'ylim',[-40 5],'ytick',-40:10:10);
set(sax(:,3),'ylim',[-100 0],'ytick',-100:25:0);
set(sax(:,4),'ylim',[150 300],'ytick',150:50:300);
set(sax(:,5),'ylim',[-0 160],'ytick',0:40:160);
for i=1:5 % align ylabels
    p=get(get(sax(1,i),'ylabel'),'position');
    set(get(sax(1,i),'ylabel'),'position',[118.1 p(2:3)])
end
orient tall;
legend({upper(name{1:2}) 'model' 'PSD ship' 'WHOI buoy'})

print -depsc2 flux_20s_oct.eps
print -dmeta flux_20s_oct.emf