-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_uvice.m
125 lines (107 loc) · 4.21 KB
/
plot_uvice.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
if ~exist('p')
p = genpath('~/MATLAB/');
addpath(p);
%load nctiles_grid in memory:
fout='~/data/geos5/MITGRID/llc90/';
grid_load(fout,5,'compact',0,0)
%displays list of grid variables:
gcmfaces_global;% disp(mygrid);
end
%V={'oceQnet','oceQsw','oceFWflx','oceSflux','oceTAUX','oceTAUY','THETA','SALT','UVEL', 'VVEL'};
V={'SIarea' 'SIheff' 'SIhsnow' 'sIceLoad' 'SIuice' 'SIvice' 'SItaux' 'SItauy'};
%var=2
%level=1
collection='iceDiag';
filo=dir(['../mit_output/' collection '*.data']);
units='[m]';
msk=mygrid.mskC(:,:,1).*mygrid.RAC;
% niter0=0;
%dt=4; %time steps per file
N=length(filo);
fignum=1;
fldTX=read_bin(['../mit_output/' filo(1).name],1,1)*0;
fldTY=read_bin(['../mit_output/' filo(1).name],1,1)*0;
area2=read_bin(['../mit_output/' filo(1).name],1,1)*0;
th=0.01;
for var=1:1
fdate=fdate0; %datetime(2000,4,14,21,0,0);
% startdate=datetime(2000,4,14,21,0,0);
% enddtate=datetime(2010,4,14,21,0,0);
i=1;
while isbetween(fdate,fdate,enddate)
if isbetween(fdate,startdate,enddate) & i<=length(filo)
area=read_bin(['../mit_output/' filo(i).name],1,1);
area(area<=th)=0;
area(area>th)=1;
area2=area2+area;
fldTX=fldTX+read_bin(['../mit_output/' filo(i).name],1,5);
fldTY=fldTY+read_bin(['../mit_output/' filo(i).name],1,6);
area2=read_bin(['../mit_output/' filo(i).name],1,1);
fdate
[V{var} ' ' num2str(i) ' ' num2str(max(fldTX/i))]
end
i=i+1;
if strcmp(DT,'hour')
fdate=fdate+hours(NDT);
elseif strcmp(DT,'day')
fdate=fdate+days(NDT);
elseif strcmp(DT,'month')
fdate=fdate+calmonths(NDT);
elseif strcmp(DT,'year')
fdate=fdate+calyears(NDT);
end
end
fldTX=fldTX/(i-1).*mygrid.mskC(:,:,1);
fldTY=fldTY/(i-1).*mygrid.mskC(:,:,1);
[fldTZ,fldTM]=calc_UEVNfromUXVY(fldTX,fldTY);
area2(area2==0)=NaN;
area2(area2>0)=1;
%zonal wind stress:
fld=fldTZ.*area2;
cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='zonal wind stress';
figureL; m_map_gcmfaces(fld,2,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
x=get(gcf,'Children');
x(2).Position=x(2).Position+[0 -0.15/2 0 0.15];
x(1).Position=x(1).Position+[0 -0.05 0 0];
y=get(gca);
y.Title.Position=y.Title.Position+[0 0.05 0];
title({'zonal ice velocity [m s^{-1}]',['mean=' num2str(round(nansum(fld.*msk)./nansum(msk)*1000)/1000) ]});
pname='figs/Ice_SIuice_np';
print(pname,'-dpng')
%meridional wind stress:
fld=fldTM.*area2;
cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/500; title0='meridional wind stress';
figureL; m_map_gcmfaces(fld,2,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
x=get(gcf,'Children');
x(2).Position=x(2).Position+[0 -0.15/2 0 0.15];
x(1).Position=x(1).Position+[0 -0.05 0 0];
y=get(gca);
y.Title.Position=y.Title.Position+[0 0.05 0];
title({'meridional ice velocity [m s^{-1}]',['mean=' num2str(round(nansum(fld.*msk)./nansum(msk)*1000)/1000) ]});
pname='figs/Ice_SIvice_np';
print(pname,'-dpng')
%zonal wind stress:
fld=fldTZ.*area2;
cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/2500; title0='zonal wind stress';
figureL; m_map_gcmfaces(fld,3,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
x=get(gcf,'Children');
x(2).Position=x(2).Position+[0 -0.15/2 0 0.15];
x(1).Position=x(1).Position+[0 -0.05 0 0];
y=get(gca);
y.Title.Position=y.Title.Position+[0 0.05 0];
title({'zonal ice velocity [m s^{-1}]',timerange,['mean=' num2str(round(nansum(fld.*msk)./nansum(msk)*1000)/1000) ]});
pname='figs/Ice_SIuice_sp';
print(pname,'-dpng')
%meridional wind stress:
fld=fldTM.*area2;
cc=[[-250:50:-100] [-75 -50] [-35:10:35] [50 75] [100:50:250]]/2500; title0='meridional wind stress';
figureL; m_map_gcmfaces(fld,3,{'myCaxis',cc},{'do_m_coast',1},{'myTitle',title0});
x=get(gcf,'Children');
x(2).Position=x(2).Position+[0 -0.15/2 0 0.15];
x(1).Position=x(1).Position+[0 -0.05 0 0];
y=get(gca);
y.Title.Position=y.Title.Position+[0 0.05 0];
title({'meridional ice velocity [m s^{-1}]',timerange,['mean=' num2str(round(nansum(fld.*msk)./nansum(msk)*1000)/1000) ]});
pname='figs/Ice_SIvice_sp';
print(pname,'-dpng')
end