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

Annotation 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 - (hide 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 gforget 1.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