/[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.4 - (hide annotations) (download)
Wed Feb 5 03:42:56 2014 UTC (11 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.3: +1 -1 lines
- rename calc_UV_div as calc_UV_conv.

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.4 seaIceDiv=1e-6*mskC.*calc_UV_conv(tmpU,tmpV); %no dh needed here
39 gforget 1.1
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 gforget 1.3 if strcmp(setDiags,'SEAICE')&myparms.diagsAreMonthly==1&userStep==-1;%computational part;
72 gforget 1.1
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 gforget 1.3 iceNorth(:,mm)=log10(mean(tmp1,2));
87 gforget 1.1 %
88     tmp1=alldiag.iceAreaSouth(:,mm:12:nt);
89     tmp1(isnan(tmp1))=0;
90 gforget 1.3 iceSouth(:,mm)=log10(mean(tmp1,2));
91 gforget 1.1 %
92     tmp1=alldiag.snowAreaNorth(:,mm:12:nt);
93     tmp1(isnan(tmp1))=0;
94 gforget 1.3 snowNorth(:,mm)=log10(mean(tmp1,2));
95 gforget 1.1 %
96     tmp1=alldiag.snowAreaSouth(:,mm:12:nt);
97     tmp1(isnan(tmp1))=0;
98 gforget 1.3 snowSouth(:,mm)=log10(mean(tmp1,2));
99 gforget 1.1 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 gforget 1.3 ' thickness distribution (in log(m$^2$))'};
123 gforget 1.1 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 gforget 1.3 ' thickness distribution (in log(m$^2$))'};
140 gforget 1.1 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 gforget 1.3 %mask out negligible amounts (e.g. due to monthly averaging)
186     msk=1+0*alldiag.seaIceConc;
187     msk(alldiag.seaIceConc<1e-2)=NaN;
188     msk(alldiag.seaIceThick<1e-2)=NaN;
189     fld=fld.*msk;
190    
191 gforget 1.1 %compute mean march and september fields
192 gforget 1.3 if seas==1;
193     tmp1=fld(:,:,3:12:240);
194     tmp1(isnan(tmp1))=0;
195     fld_seas=mean(tmp1,3);
196     else;
197     tmp1=fld(:,:,9:12:240);
198     tmp1(isnan(tmp1))=0;
199     fld_seas=mean(tmp1,3);
200 gforget 1.1 end;
201     fld_seas(fld_seas==0)=NaN;
202    
203     %ice concentration
204     figureL; set(gcf,'Renderer','zbuffer');
205     cc=list_cc{vv}; if doAnomalies; cc=list_ccAno{vv}; end;
206     caxi={'myCaxis',cc}; docb={'doCbar',1};
207     m_map_gcmfaces(fld_seas,pp,caxi,docb);
208     myCaption={myYmeanTxt,mon,' mean -- ',list_tit{vv}};
209     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
210    
211     end;%for vv=1:length(list_var);
212     end;%for seas=1:2;
213     end;%for hemi=1:2;
214    
215     end;
216    

  ViewVC Help
Powered by ViewVC 1.1.22