/[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.2 - (hide annotations) (download)
Fri Jan 31 23:21:35 2014 UTC (11 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.1: +15 -9 lines
- add subsections
- fix 'divergence'

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 gforget 1.2 %convergence in kiloton/s
38 gforget 1.1 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 gforget 1.2 if addToTex; write2tex(fileTex,1,'Monthly Thickness Distribution',2); end;
76    
77 gforget 1.1 %compute seasonal cycle in thickness distributions
78    
79     y=alldiag.iceThickNorth(:,1:12);
80     x=ones(size(y,1),1)*[1:12];
81     iceNorth=zeros(size(x)); iceSouth=zeros(size(x));
82     snowNorth=zeros(size(x)); snowSouth=zeros(size(x));
83     for mm=1:12;
84     tmp1=alldiag.iceAreaNorth(:,mm:12:nt);
85     tmp1(isnan(tmp1))=0;
86     iceNorth(:,mm)=log10(median(tmp1,2));
87     %
88     tmp1=alldiag.iceAreaSouth(:,mm:12:nt);
89     tmp1(isnan(tmp1))=0;
90     iceSouth(:,mm)=log10(median(tmp1,2));
91     %
92     tmp1=alldiag.snowAreaNorth(:,mm:12:nt);
93     tmp1(isnan(tmp1))=0;
94     snowNorth(:,mm)=log10(median(tmp1,2));
95     %
96     tmp1=alldiag.snowAreaSouth(:,mm:12:nt);
97     tmp1(isnan(tmp1))=0;
98     snowSouth(:,mm)=log10(median(tmp1,2));
99     end;
100     iceNorth(~isfinite(iceNorth))=NaN;
101     iceSouth(~isfinite(iceSouth))=NaN;
102     snowNorth(~isfinite(snowNorth))=NaN;
103     snowSouth(~isfinite(snowSouth))=NaN;
104    
105     %now display seasonal cycle
106     x=[x x+12]-0.5; y=[y y]; z=[iceNorth iceNorth];
107    
108     figureL; set(gcf,'Renderer','zbuffer');
109     subplot(2,1,1);
110     z=[iceNorth iceNorth]; pcolor(x,y,z); shading interp;
111     axis([1 13 0 6]); grid on; caxis([9.5 12]); colormap(jet(25));
112     colorbar; xlabel('month of year'); ylabel('log(thickness)');
113     title('log(ice area) -- North. Hemi.');
114     subplot(2,1,2);
115     z=[snowNorth snowNorth]; pcolor(x,y,z); shading interp;
116     axis([1 13 0 1.2]); grid on; caxis([9.5 12]+0.5); colormap(jet(25));
117     colorbar; xlabel('month of year'); ylabel('log(thickness)');
118     title('log(snow area) -- North. Hemi.');
119    
120     myCaption={myYmeanTxt,'Northern Hemisphere :',...
121     'monthly mean ice (top) and snow (bottom)',...
122     ' thickness distribution (in log(m2))'};
123     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
124    
125     figureL; set(gcf,'Renderer','zbuffer');
126     subplot(2,1,1);
127     z=[iceSouth iceSouth]; pcolor(x,y,z); shading interp;
128     axis([1 13 0 4]); grid on; caxis([9.5 12]); colormap(jet(25));
129     colorbar; xlabel('month of year'); ylabel('log(thickness)');
130     title('log(ice area) -- South. Hemi.');
131     subplot(2,1,2);
132     z=[snowSouth snowSouth]; pcolor(x,y,z); shading interp;
133     axis([1 13 0 1.2]); grid on; caxis([9.5 12]+0.5); colormap(jet(30));
134     colorbar; xlabel('month of year'); ylabel('log(thickness)');
135     title('log(snow area) -- South. Hemi.');
136    
137     myCaption={myYmeanTxt,'Southern Hemisphere :',...
138     'monthly mean ice (top) and snow (bottom)',...
139     ' thickness distribution (in log(m2))'};
140     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
141    
142     end;
143    
144     %maps of mean march and september ice fields
145    
146     list_var={'seaIceConc',...
147     'seaIceThick',...
148     'seaSnowThick',...
149     'seaIceStrFun',...
150     'seaIceDiv'};
151    
152     list_cc={[0:0.1:1],...
153     [0:0.1:0.5 0.75 1:0.5:2.5 3:5],...
154     [0:0.1:0.5 0.75 1:0.5:2.5 3:5]/5,...
155     1/50*[-5:-3 -2.5:0.5:-1 -0.5:0.2:0.5 1:0.5:2.5 3:5],...
156     1/5*[-5:-3 -2.5:0.5:-1 -0.5:0.2:0.5 1:0.5:2.5 3:5]};
157    
158     list_ccAno={2*[-0.2:0.04:0.2],...
159     4*[-0.2:0.04:0.2],...
160     [-0.2:0.04:0.2],...
161     0.25*[-0.2:0.04:0.2],...
162     [-0.2:0.04:0.2]};
163    
164     list_tit={' ice concentration (unitless)',...
165     ' ice thickness (m)',...
166     ' snow thickness (m)',...
167     ' ice+snow streamfunction (megaton/s)',...
168 gforget 1.2 ' ice+snow convergence (kiloton/s)'};
169 gforget 1.1
170     for hemi=1:2;
171     for seas=1:2;
172 gforget 1.2
173     %select month
174     if seas==1; mon='March'; else; mon='September'; end;
175    
176     %select projection
177     if hemi==1; pp=2; txt='Northern '; else; pp=3; txt='Southern '; end;
178    
179     if addToTex; write2tex(fileTex,1,[txt 'Hem. in ' mon],2); end;
180    
181 gforget 1.1 for vv=1:length(list_var);
182    
183     eval(['fld=alldiag.' list_var{vv} ';']);
184    
185     %compute mean march and september fields
186 gforget 1.2 if seas==1; fld_seas=nanmedian(fld(:,:,3:12:240),3);
187     else; fld_seas=nanmedian(fld(:,:,9:12:240),3);
188 gforget 1.1 end;
189     fld_seas(fld_seas==0)=NaN;
190    
191     %ice concentration
192     figureL; set(gcf,'Renderer','zbuffer');
193     cc=list_cc{vv}; if doAnomalies; cc=list_ccAno{vv}; end;
194     caxi={'myCaxis',cc}; docb={'doCbar',1};
195     m_map_gcmfaces(fld_seas,pp,caxi,docb);
196     myCaption={myYmeanTxt,mon,' mean -- ',list_tit{vv}};
197     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
198    
199     end;%for vv=1:length(list_var);
200     end;%for seas=1:2;
201     end;%for hemi=1:2;
202    
203     end;
204    

  ViewVC Help
Powered by ViewVC 1.1.22