Repository navigation
Expand file tree
/
Copy pathmake_FV_3Layers_subset.m
More file actions
169 lines (153 loc) · 5.6 KB
/
Copy pathmake_FV_3Layers_subset.m
File metadata and controls
169 lines (153 loc) · 5.6 KB
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
%function[]=make_FV_3Layers_subset();
%(c) dmitry.aleynik@sams.ac.uk 2017.09.01 &
% __o_O__¬ FASTNEt &
% /_____/ ~~~~~~<@})))< ~ ~ ~~~~~ ~~~ SAMS Glider mission 4, Malin Shelf &
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
if 1,...
clear all
close all
this_dir=pwd; DRV=this_dir(1:1);
Tb=datenum('2017-09-06 08:00:00',31); % dye release 1kg rodamin
Te=datenum('2017-09-07 16:00:00',31);
tic;
addpath M:/Mar_Phys/matlab/m_map
addpath M:/Mar_Phys/matlab/seawater
addpath M:/Mar_Phys/matlab/general/
a=load('etive6_0coast.mat');
path_fig='../fig/';
path_mat='../mat/';
filePO=['Poltips_Oban_20170901-0930.txt'];
[PO] = import_pol_oban(filePO);
v_po = datevec(PO.mtime(1) );
time_offset = 678942; % from FVCOM time to matlab time
[IM_range_lat ,IM_range_lon] = deal([56.445 56.496667 ],[ -(5+28.0/60) -(5+22.5/60)]); cm='_AM'; %Ardmucknish
ml2m=1852*60; rlat= mean(IM_range_lat);
xv=[IM_range_lon]; yv=[IM_range_lat];
m_proj('UTM','lon',[ IM_range_lon],'lat',[IM_range_lat],'zon',30,'ell','wgs84')
load('A_Grid.mat','-mat');
%% ======================
%% Get gridded subsets
now1=now; disp([ ' Loading Start at: ' datestr(now1,31)]);
dir_fvmg=dir ([path_mat 'FV_Grd20*.mat']);
% vars={'el';'t';'s'; 'u';'v';'ww';'ua';'va';'uwind_speed';'vwind_speed'};
vars={'el';'u';'v';'ww';};
[fn, fi]=size(dir_fvmg);
% fn=10; % 30 h, the latest files 11-13 are slow to get at laptop
% fn=13; % 32 h, 48.230603 seconds at sams PC
for fi=1:fn,...
A=[];got=[];
A=load([path_mat dir_fvmg(fi).name ]);
got=[]; got= find(A.FV.mt >=Tb & A.FV.mt <=Te);
if isempty(got), continue;end
%% define vertical levels etc :
if fi==1,...
FV=[]; mt=[];
FV.mask=A.FV.mask; FV.x=A.FV.x; FV.y=A.FV.y; FV.siglev=A.FV.siglev;
FV.h=A.FV.h; mt=A.FV.mt;
% [ nx,ny] = size(A.FV.h);
% [nt, nx,ny] = size(A.FV.el);
[nt,nz,nx,ny] = size(A.FV.u);
nt = length(got);
FV.mt=A.FV.mt;
FV.el=A.FV.el;
FV.u=A.FV.u;
FV.v=A.FV.v;
FV.ww=A.FV.ww;
FV.fi(fi,1)=fi;
FV.fnm(fi,:)=dir_fvmg(fi).name ;
FV.go_nt(fi,1)=length(got);
else
go=[]; go = find(A.FV.mt > FV.mt(end));
FV.mt=cat(1,FV.mt,A.FV.mt(go) );
FV.el=cat(1,FV.el,A.FV.el(go,:,:) );
FV.u =cat(1,FV.u ,A.FV.u( go,:,:,:));
FV.v =cat(1,FV.v ,A.FV.v( go,:,:,:));
FV.ww=cat(1,FV.ww,A.FV.ww(go,:,:,:));
FV.fi(fi,1)=fi;
FV.go_nt(fi,1)=length(go);
FV.fnm(fi,:)=dir_fvmg(fi).name ;
end
disp([ ' Loading ' num2str(fi) ' ' dir_fvmg(fi).name ' at: ' datestr(now,31),' dt=' num2str((now-now1)*24*3600) ' s']);
end; %fn
toc ; % ~100 seconds over 10 files , 30 muntes for 13 files
disp([ ' Loading End at: ' datestr(now,31),' dt=' num2str((now-now1)*24*3600) ' sec']);
clear A;
%% get the vertical sigma-layers depth : zs
nt=length(FV.mt); % tic
[nz ~]=size(FV.siglev); nz=nz-1;
yz=[]; zs=[]; zs=zeros(nt,nz+1,nx,ny);
for it=1:nt,...
for i=1:nx,
for j=1:ny,...
if FV.mask(i,j)==0, continue; end
% yz=zeros(nz+1,1);
yz=-( -( FV.h(i,j)+FV.el(it,i,j) ).*FV.siglev(1:end ) ) + FV.el(it,i,j);
% -( -( 36 +2 ) * [0:-1] ) + 2
zs(it,:,i,j) =yz;
end;
end;
end
FV.zs=zs;
clear zs; % toc % 56 sec;
disp([ ' Zlayers End at: ' datestr(now,31),' dt=' num2str((now-now1)*24*3600) ' sec']);
end;
FV3=FV;
FV3=rmfield(FV3,{'u';'v';'ww';'zs'});
FV3.layers=[1,fix(nz/2),nz];
for iz=1:length(FV3.layers),...
jz=FV3.layers(iz);
% a=[]; a=FV.zs(:,jz,:,:);
FV3.zl(:,iz,:,:)=FV.zs(:,jz,:,:); ... end;
FV3.u( :,iz,:,:)=FV.u(:,jz,:,:);
FV3.v( :,iz,:,:)=FV.v(:,jz,:,:);
FV3.ww( :,iz,:,:)=FV.ww(:,jz,:,:);
end
dir_fvmg=dir ([path_mat 'FV_Grd20*.mat']);
db=datestr(FV3.mt(1),'yyyymmdd-HHMM');
de=datestr(FV3.mt(end),'yyyymmdd-HHMM');
save([path_mat 'FV_3Layers_Grd_' db '-' de '.mat'],'FV3','-v7.3');
disp([ ' 3 Layers saved at: ' datestr(now,31),' dt=' num2str((now-now1)*24*3600) ' sec']);
% test it here
a=FV3.zl(:,:,10,10);
b=FV.zs(:,:,10,10);
plot(FV3.mt,b); hold on;
plot(FV3.mt,a,'linewidth',2); hold on;
plot(FV3.mt,FV3.el(:,10,10),'k'); hold on;
set(gca,'xlim',[min(FV.mt) max(FV.mt)]);
ylabel('depth,m');
set(gca,'YMinorTick','on','XMinorTick','on');
datetick('x','keeplimits');
print(['-f'],'-dpng','-loose','-r300',['z_sigma_i10j10' '.png']);
if 0
%% remake FV-3hours files got=[]; got= find(FV.mt >=Tb & FV.mt <=Te);
nt=length(got);mt=FV.mt;
got0 = find(mod(mt(:),3/24)==0);
got0(end+1)=got(end);
nk=length(got0);
FV0=FV; vars={'el';'u';'v';'ww';};
for ik=1:nk,...
if ik==1,
ib=1;
else
ib=got0(ik-1)+1;
end
ie=got0(ik);
go=[ib:ie]'; ng =length(go);
FV=[];
FV=FV0;
FV=rmfield(FV,{'u';'v';'ww';'zs';'el',;'mt'});
FV.mt(1:ng,1) =FV0.mt(go,1) ;
FV.el(1:ng,:,: )=FV0.el(go,:,: );
FV.u( 1:ng,:,:,:)=FV0.u( go,:,:,:);
FV.v( 1:ng,:,:,:)=FV0.v( go,:,:,:);
FV.ww(1:ng,:,:,:)=FV0.ww(go,:,:,:);
% FV.zs(1:ng,:,:,:)=FV0.zs(go,:,:,:);
db=datestr(FV.mt(1),'yyyymmdd-HHMM');
de=datestr(FV.mt(end),'yyyymmdd-HHMM');
FVgrd_nm=['FV_Grd' db '_' de];
save([path_mat FVgrd_nm '.mat'],'FV','-v7.3');
disp(['saved :' FVgrd_nm ', nt=', num2str(length(FV.mt)),'; ' , datestr(now,31)]);
end
end
% plt_FV_subset
% plt_FV_subset_vect