/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_SEAICE.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_SEAICE.m

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


Revision 1.1 - (hide annotations) (download)
Tue Jan 28 17:07:13 2014 UTC (11 years, 6 months ago) by gforget
Branch: MAIN
- diags_display.m : omit annual running mean when diagsNbRec<=12
- diags_pre_process.m and diags_set_E.m : fix typos
- diags_set_D.m and diags_set_E.m : alow for 'mat/' subdirectories
- added : diags_set_SEAICE.m dedicated to seaice and snow diagnostics

1 gforget 1.1
2     if strcmp(setDiags,'SEAICE')&userStep==1;%diags to be computed
3     listDiags='seaIceConc seaIceThick seaSnowThick seaIceStrFun seaIceDiv';
4     listDiags=[listDiags ' iceThickNorth iceThickSouth'];
5     listDiags=[listDiags ' iceAreaNorth iceAreaSouth'];
6     listDiags=[listDiags ' snowThickNorth snowThickSouth'];
7     listDiags=[listDiags ' snowAreaNorth snowAreaSouth'];
8     end;
9    
10     if strcmp(setDiags,'SEAICE')&userStep==2;%input files and variables
11     listFlds={'SIarea ','SIheff ','SIhsnow ',...
12     'ADVxHEFF','ADVyHEFF','ADVxSNOW','ADVySNOW',...
13     'DFxEHEFF','DFyEHEFF','DFxESNOW','DFyESNOW'};
14     listFldsNames=deblank(listFlds);
15     listFiles={'state_2d_set1'};
16     listSubdirs={[dirModel 'diags/STATE/' ]};
17     end;
18    
19     if strcmp(setDiags,'SEAICE')&userStep==3;%computational part;
20    
21     %ice masks
22     mskC=mygrid.mskC(:,:,1); mskC(SIheff==0)=NaN;
23     mskW=mygrid.mskW(:,:,1); mskW(ADVxHEFF==0)=NaN;
24     mskS=mygrid.mskS(:,:,1); mskS(ADVyHEFF==0)=NaN;
25    
26     %ice and snow thickness
27     seaIceConc=mskC.*SIarea;
28     seaIceThick=SIheff./seaIceConc;
29     seaSnowThick=SIhsnow./seaIceConc;
30    
31     %transports in kg/s
32     tmpU=(myparms.rhoi*DFxEHEFF+myparms.rhosn*DFxESNOW+...
33     myparms.rhoi*ADVxHEFF+myparms.rhosn*ADVxSNOW);
34     tmpV=(myparms.rhoi*DFyEHEFF+myparms.rhosn*DFyESNOW+...
35     myparms.rhoi*ADVyHEFF+myparms.rhosn*ADVySNOW);
36    
37     %divergence in kiloton/s
38     seaIceDiv=1e-6*mskC.*calc_UV_div(tmpU,tmpV); %no dh needed here
39    
40     %streamfunction in kg/s
41     %(factor 1e6 to compensare for calc_barostream internal conversion)
42     seaIceStrFun=1e6*calc_barostream(tmpU.*mskW,tmpV.*mskS,1,{});
43     %convert to megaton/s
44     seaIceStrFun=1e-9*mskC.*seaIceStrFun;
45    
46     %thickness distributions per hemisphere
47     xx=convert2vector(seaIceThick);
48     yy=convert2vector(1+0*seaIceThick);
49     zz=convert2vector(seaIceConc.*mygrid.RAC);
50    
51     mm=convert2vector(mygrid.YC>0); jj=find(mm.*zz>0);
52     [x,y,z,n]=MITprof_stats(xx(jj),[0:0.05:10],yy(jj),[0 1 2],'sum',zz(jj));
53     iceThickNorth=x(:,2); iceAreaNorth=z(:,2);
54    
55     mm=convert2vector(mygrid.YC<0); jj=find(mm.*zz>0);
56     [x,y,z,n]=MITprof_stats(xx(jj),[0:0.05:10],yy(jj),[0 1 2],'sum',zz(jj));
57     iceThickSouth=x(:,2); iceAreaSouth=z(:,2);
58    
59     xx=convert2vector(seaSnowThick);
60    
61     mm=convert2vector(mygrid.YC>0); jj=find(mm.*zz>0);
62     [x,y,z,n]=MITprof_stats(xx(jj),[0:0.05:10],yy(jj),[0 1 2],'sum',zz(jj));
63     snowThickNorth=x(:,2); snowAreaNorth=z(:,2);
64    
65     mm=convert2vector(mygrid.YC<0); jj=find(mm.*zz>0);
66     [x,y,z,n]=MITprof_stats(xx(jj),[0:0.05:10],yy(jj),[0 1 2],'sum',zz(jj));
67     snowThickSouth=x(:,2); snowAreaSouth=z(:,2);
68    
69     end;
70    
71     if strcmp(setDiags,'SEAICE')&userStep==-1;%computational part;
72    
73     if ~doAnomalies;
74    
75     %compute seasonal cycle in thickness distributions
76    
77     y=alldiag.iceThickNorth(:,1:12);
78     x=ones(size(y,1),1)*[1:12];
79     iceNorth=zeros(size(x)); iceSouth=zeros(size(x));
80     snowNorth=zeros(size(x)); snowSouth=zeros(size(x));
81     for mm=1:12;
82     tmp1=alldiag.iceAreaNorth(:,mm:12:nt);
83     tmp1(isnan(tmp1))=0;
84     iceNorth(:,mm)=log10(median(tmp1,2));
85     %
86     tmp1=alldiag.iceAreaSouth(:,mm:12:nt);
87     tmp1(isnan(tmp1))=0;
88     iceSouth(:,mm)=log10(median(tmp1,2));
89     %
90     tmp1=alldiag.snowAreaNorth(:,mm:12:nt);
91     tmp1(isnan(tmp1))=0;
92     snowNorth(:,mm)=log10(median(tmp1,2));
93     %
94     tmp1=alldiag.snowAreaSouth(:,mm:12:nt);
95     tmp1(isnan(tmp1))=0;
96     snowSouth(:,mm)=log10(median(tmp1,2));
97     end;
98     iceNorth(~isfinite(iceNorth))=NaN;
99     iceSouth(~isfinite(iceSouth))=NaN;
100     snowNorth(~isfinite(snowNorth))=NaN;
101     snowSouth(~isfinite(snowSouth))=NaN;
102    
103     %now display seasonal cycle
104     x=[x x+12]-0.5; y=[y y]; z=[iceNorth iceNorth];
105    
106     figureL; set(gcf,'Renderer','zbuffer');
107     subplot(2,1,1);
108     z=[iceNorth iceNorth]; pcolor(x,y,z); shading interp;
109     axis([1 13 0 6]); grid on; caxis([9.5 12]); colormap(jet(25));
110     colorbar; xlabel('month of year'); ylabel('log(thickness)');
111     title('log(ice area) -- North. Hemi.');
112     subplot(2,1,2);
113     z=[snowNorth snowNorth]; pcolor(x,y,z); shading interp;
114     axis([1 13 0 1.2]); grid on; caxis([9.5 12]+0.5); colormap(jet(25));
115     colorbar; xlabel('month of year'); ylabel('log(thickness)');
116     title('log(snow area) -- North. Hemi.');
117    
118     myCaption={myYmeanTxt,'Northern Hemisphere :',...
119     'monthly mean ice (top) and snow (bottom)',...
120     ' thickness distribution (in log(m2))'};
121     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
122    
123     figureL; set(gcf,'Renderer','zbuffer');
124     subplot(2,1,1);
125     z=[iceSouth iceSouth]; pcolor(x,y,z); shading interp;
126     axis([1 13 0 4]); grid on; caxis([9.5 12]); colormap(jet(25));
127     colorbar; xlabel('month of year'); ylabel('log(thickness)');
128     title('log(ice area) -- South. Hemi.');
129     subplot(2,1,2);
130     z=[snowSouth snowSouth]; pcolor(x,y,z); shading interp;
131     axis([1 13 0 1.2]); grid on; caxis([9.5 12]+0.5); colormap(jet(30));
132     colorbar; xlabel('month of year'); ylabel('log(thickness)');
133     title('log(snow area) -- South. Hemi.');
134    
135     myCaption={myYmeanTxt,'Southern Hemisphere :',...
136     'monthly mean ice (top) and snow (bottom)',...
137     ' thickness distribution (in log(m2))'};
138     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
139    
140     end;
141    
142     %maps of mean march and september ice fields
143    
144     list_var={'seaIceConc',...
145     'seaIceThick',...
146     'seaSnowThick',...
147     'seaIceStrFun',...
148     'seaIceDiv'};
149    
150     list_cc={[0:0.1:1],...
151     [0:0.1:0.5 0.75 1:0.5:2.5 3:5],...
152     [0:0.1:0.5 0.75 1:0.5:2.5 3:5]/5,...
153     1/50*[-5:-3 -2.5:0.5:-1 -0.5:0.2:0.5 1:0.5:2.5 3:5],...
154     1/5*[-5:-3 -2.5:0.5:-1 -0.5:0.2:0.5 1:0.5:2.5 3:5]};
155    
156     list_ccAno={2*[-0.2:0.04:0.2],...
157     4*[-0.2:0.04:0.2],...
158     [-0.2:0.04:0.2],...
159     0.25*[-0.2:0.04:0.2],...
160     [-0.2:0.04:0.2]};
161    
162     list_tit={' ice concentration (unitless)',...
163     ' ice thickness (m)',...
164     ' snow thickness (m)',...
165     ' ice+snow streamfunction (megaton/s)',...
166     ' ice+snow divergence (kiloton/s)'};
167    
168     for hemi=1:2;
169     for seas=1:2;
170     for vv=1:length(list_var);
171    
172     eval(['fld=alldiag.' list_var{vv} ';']);
173    
174     %compute mean march and september fields
175     if seas==1; fld_seas=nanmedian(fld(:,:,3:12:240),3); mon='March';
176     else; fld_seas=nanmedian(fld(:,:,9:12:240),3); mon='September';
177     end;
178     fld_seas(fld_seas==0)=NaN;
179    
180     %choose projection
181     if hemi==1; pp=2; txt='Northern ';
182     else; pp=3; txt='Southern ';
183     end;
184    
185     %ice concentration
186     figureL; set(gcf,'Renderer','zbuffer');
187     cc=list_cc{vv}; if doAnomalies; cc=list_ccAno{vv}; end;
188     caxi={'myCaxis',cc}; docb={'doCbar',1};
189     m_map_gcmfaces(fld_seas,pp,caxi,docb);
190     myCaption={myYmeanTxt,mon,' mean -- ',list_tit{vv}};
191     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
192    
193     end;%for vv=1:length(list_var);
194     end;%for seas=1:2;
195     end;%for hemi=1:2;
196    
197     end;
198    

  ViewVC Help
Powered by ViewVC 1.1.22