1 |
function []=bp_ECCO_GRACE(); |
2 |
% BP_ECCO_GRACE |
3 |
% - reads ECCO and GRACE bottom pressure fields from file |
4 |
% - computes monthly area average over Barents Sea |
5 |
% - display monthly time series and correlation |
6 |
|
7 |
fprintf('\n'); help bp_ECCO_GRACE; |
8 |
|
9 |
%load grid from nctiles_grid/ |
10 |
grid_load; gcmfaces_global; |
11 |
|
12 |
%option to subtract ECCO time-variable global-mean BP (by default rmGloMean=0 |
13 |
% and GRACE_jpl_rl05m_SpatialMean.asc is instead added to GRACE BP maps) |
14 |
rmGloMean=0; |
15 |
|
16 |
%input directories |
17 |
dir_GRACE='input_ecco/input_bp/'; |
18 |
dir_ECCO='nctiles_monthly/PHIBOT/'; |
19 |
|
20 |
%read GRACE global mean |
21 |
glo_GRACE=NaN*zeros(1,240); |
22 |
fid=fopen([dir_GRACE 'GRACE_jpl_rl05m_SpatialMean.asc'],'rt'); |
23 |
for ii=1:4; tmp1=fgetl(fid); end; |
24 |
for ii=1:240; |
25 |
tmp1=fgetl(fid); |
26 |
tmp1=str2num(tmp1); |
27 |
glo_GRACE(ii)=tmp1(4); |
28 |
end; |
29 |
fclose(fid); |
30 |
glo_GRACE(glo_GRACE==-999)=NaN; |
31 |
|
32 |
%read GRACE |
33 |
fld_GRACE=repmat(NaN*mygrid.XC,[1 1 240]); |
34 |
for yy=1992:2011; |
35 |
tmp1=read_bin([dir_GRACE 'GRACE_jpl_rl05m_' num2str(yy)]); |
36 |
fld_GRACE(:,:,[1:12]+12*(yy-1992))=tmp1; |
37 |
end; |
38 |
|
39 |
%read ECCO |
40 |
fld_ECCO=read_nctiles([dir_ECCO 'PHIBOT']); |
41 |
|
42 |
%convert ECCO to cm |
43 |
fld_ECCO=10.1937*fld_ECCO; |
44 |
|
45 |
%apply common NaN mask to ECCO and GRACE maps |
46 |
msk=1*(fld_GRACE>-999); msk(msk==0)=NaN; |
47 |
fld_GRACE=msk.*fld_GRACE; |
48 |
fld_ECCO=msk.*fld_ECCO; |
49 |
|
50 |
%NaN mask for GRACE / ECCO time series |
51 |
msk_tim=1*(nansum(msk,0)>0); |
52 |
msk_tim(msk_tim==0)=NaN; |
53 |
|
54 |
%subtract time mean |
55 |
fld_GRACE=fld_GRACE-repmat(nanmean(fld_GRACE,3),[1 1 240]); |
56 |
fld_ECCO=fld_ECCO-repmat(nanmean(fld_ECCO,3),[1 1 240]); |
57 |
|
58 |
%add global mean for GRACE |
59 |
fld_GRACE=fld_GRACE+mk3D(glo_GRACE,fld_GRACE); |
60 |
|
61 |
%subtract time variable global mean |
62 |
if rmGloMean; |
63 |
tmp1=nanmedian(fld_GRACE,[],0); |
64 |
fld_GRACE=fld_GRACE-mk3D(tmp1,fld_GRACE); |
65 |
tmp1=nanmedian(fld_ECCO,[],0); |
66 |
fld_ECCO=fld_ECCO-mk3D(tmp1,fld_ECCO); |
67 |
end; |
68 |
|
69 |
%Barents Sea (area weighted mean) time series: |
70 |
msk=v4_basin('barents').*mygrid.mskC(:,:,1);; |
71 |
mskXrac=repmat(msk.*mygrid.RAC,[1 1 240]); |
72 |
bp_ECCO=nansum(fld_ECCO.*mskXrac,0)./nansum(mskXrac,0); |
73 |
bp_GRACE=nansum(fld_GRACE.*mskXrac,0)./nansum(mskXrac,0); |
74 |
|
75 |
%apply common NaN mask to ECCO and GRACE time series |
76 |
bp_GRACE=msk_tim.*bp_GRACE; |
77 |
bp_ECCO=msk_tim.*bp_ECCO; |
78 |
|
79 |
%display results |
80 |
figure; set(gca,'FontSize',12); x_tim=[0.5:239.5]/12+1992; |
81 |
plot(x_tim,bp_GRACE,'b.-'); hold on; plot(x_tim,bp_ECCO,'r.-'); |
82 |
legend('GRACE','ECCO'); ylabel('Bottom pressure anomaly (in cm)'); |
83 |
title('Barents Sea, area averaged, monthly means'); |
84 |
|
85 |
tmp0=~isnan(msk_tim); |
86 |
tmp1=corrcoef(bp_ECCO(tmp0),bp_GRACE(tmp0)); |
87 |
fprintf(' Correlation coefficient = %g\n\n',tmp1(2,1)); |
88 |
|