start='June';
lead='L/0.5/3.5/';

%or

start='May';
lead='L/01.5/4.5/';

loaddap(['http://awr:DLaccess@iridl.ldeo.columbia.edu/expert/home/.mbell/SOURCES/.IMD/.NCC1-2005/.v2p0/.rf/T/(Jun-Sep)seasonalAverage/T/(Jun-Sep%201982)(Jun-Sep%202004)RANGE/T//S/renameGRID/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.MONTHLY/.prec/S/(0000%201%20',start,'%201982-2009)VALUES/',lead,'RANGE%5BM/L%5Daverage%5BX/Y%5DregridLinear/X/68/90/RANGE/Y/5/35/RANGE/S/3/shiftGRID/8640000/mul/dods'])

ECHAM4p5_CA=prec.prec;

loaddap(['http://awr:DLaccess@iridl.ldeo.columbia.edu/expert/home/.mbell/SOURCES/.IMD/.NCC1-2005/.v2p0/.rf/T/(Jun-Sep)seasonalAverage/T/(Jun-Sep%201982)(Jun-Sep%202004)RANGE/T//S/renameGRID/SOURCES/.IRI/.FD/.ECHAM4p5/.Forecast/.ca_sst/.ensemble24/.NEWMONTHLY/.may2009/.surface/.prcp/T/%28Jun%202009%29%28Sep%202009%29RANGEEDGES/%5BM/T%5Daverage%5BX/Y%5DregridLinear/X/68/90/RANGE/Y/5/35/RANGE/S/3/shiftGRID/8640000/mul/dods'])

ECHAM4p5_CA(:,:,28) = prcp.prcp;


loaddap(['http://awr:DLaccess@iridl.ldeo.columbia.edu/expert/home/.mbell/SOURCES/.IMD/.NCC1-2005/.v2p0/.rf/T/(Jun-Sep)seasonalAverage/T/(Jun-Sep%201982)(Jun-Sep%202004)RANGE/T//S/renameGRID/SOURCES/.IRI/.MP/.RESEARCH/.COUPLED/.GLOBAL/.ECHAM4p5-GML/.CFS_sst/.ensemble12/.MONTHLY/.Surface/.prec/S/(0000%201%20',start,'%201982-2009)VALUES/',lead,'RANGE%5BM/L%5Daverage%5BX/Y%5DregridLinear/X/68/90/RANGE/Y/5/35/RANGE/S/3/shiftGRID/8640000/mul/dods'])

ECHAM4p5_GML=prec.prec;

loaddap(['http://awr:DLaccess@iridl.ldeo.columbia.edu/expert/home/.mbell/SOURCES/.IMD/.NCC1-2005/.v2p0/.rf/T/(Jun-Sep)seasonalAverage/T/(Jun-Sep%201982)(Jun-Sep%202004)RANGE/T//S/renameGRID/SOURCES/.NOAA/.NCEP/.EMC/.CFS/.MONTHLY/.surface/.prate/S/(0000%201%20',start,'%201982-2009)VALUES/',lead,'RANGE%5BM/L%5Daverage%5BX/Y%5DregridLinear/X/68/90/RANGE/Y/5/35/RANGE/S/3/shiftGRID/8640/mul/dods'])
CFS = prate.prate;

loaddap('http://awr:DLaccess@iridl.ldeo.columbia.edu/expert/home/.mbell/SOURCES/.IMD/.NCC1-2005/.v2p0/.rf/T/(Jun-Sep)seasonalAverage/T/(Jun-Sep%201982)(Jun-Sep%202004)RANGE/X/67.5/90.5/RANGEEDGES/Y/6.5/35.5/RANGEEDGES/dods')
OBS=rf.rf;

nobs = 23;
nf = 28;
  
mm = ECHAM4p5_GML+ECHAM4p5_CA;
mmm = mean(mm(:,:,1:nobs),3);
mma = mm - mmm(:,:,ones(1,nf));
mm_std = std(mma(:,:,1:nobs),[],3);
mma = mma./mm_std(:,:,ones(1,nf));
c=correlation2D(mm(:,:,1:nobs),OBS);

subplot(2,2,1) % #1 of a 2 x 2 set of plots
my_pcolor(X,Y,c) %color plot
cbar('vert',[-1:0.25:1]) %color bar
pacific_noresize_at_all%coast line     
title(['GML correlation'])

r=c;
r(r<0.1)=eps;;
sf = 1./sqrt(r.^2./(1-r.^2));
%
b13=norminv(1/3)*sqrt(sf.^2+1);
b23=norminv(2/3)*sqrt(sf.^2+1);
Pb = normcdf(b13,mma(:,:,end),sf);
Pa = 1-normcdf(b23,mma(:,:,end),sf);
Pn = 1-Pa-Pb;


subplot(2,2,2)
my_pcolor(X,Y,100*Pb)
cbar('vert',[15,20,25,30,35,40,45,50,55]-2)
pacific_noresize_at_all     
title('P(below)')

subplot(2,2,3)
my_pcolor(X,Y,100*Pn)
cbar('vert',[15,20,25,30,35,40,45,50,55]-2)
pacific_noresize_at_all     
title('P(normal)')

subplot(2,2,4)
my_pcolor(X,Y,100*Pa)
cbar('vert',[15,20,25,30,35,40,45,50,55]-2)
pacific_noresize_at_all     
title('P(above)')
