/[MITgcm]/MITgcm_contrib/gael/matlab_class/ecco_v4/bp_ECCO_GRACE.m
ViewVC logotype

Contents of /MITgcm_contrib/gael/matlab_class/ecco_v4/bp_ECCO_GRACE.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Tue Apr 4 15:41:35 2017 UTC (8 years, 3 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66f, HEAD
- Example comparison between ECCO and GRACE bottom pressure time series

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

  ViewVC Help
Powered by ViewVC 1.1.22