-
Notifications
You must be signed in to change notification settings - Fork 2
/
calculate_G_v_size.m
123 lines (121 loc) · 3.09 KB
/
calculate_G_v_size.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
% calculate_G_v_size - revised version of calculate_G that assumes
% plot_floc variables exist and uses native Matlab
load('ustar_av')
% url=['http://geoport.whoi.edu/thredds/dodsC/clay/usgs/users/aretxabaleta/MVCO/ocean_his_',cas,'.nc'];
% nc=ncgeodataset(url);
% %
% tim=nc{'ocean_time'}(:)/3600/24+ustar_av.dn(1);
% nt=length(tim);
% zeta = squeeze( nc{'zeta'}(:,3,3) );
% h = nc{'h'}(3,3);
% Cs_r = nc{'Cs_r'}(:) ; Cs_w = nc{'Cs_w'}(:) ;
% nz=length(Cs_r);
% s_w = nc{'s_w'}(:) ; s_rho = nc{'s_rho'}(:) ; hc = nc{'hc'}(:) ;
% for k = 1:nz+1
% zw(:,k) = zeta.*(1+s_w(k)*hc./h-hc*Cs_w(k)./h+Cs_w(k))+(h - hc)*Cs_w(k);
% end
% for k = 1:nz
% z_r(:,k) = zeta.*(1+s_rho(k)*hc./h-hc*Cs_r(k)./h+Cs_r(k))+hc*s_rho(k)+(h - hc)*Cs_r(k);
% end
% dz = diff(zw,1,2);
i=3; j=4;
% rho=nanmean(nanmean(nc{'rho'}(:,:,3,3)))+1000;
% bstru=nc{'bustrcwmax'}(:,3,3);
% bstrv=nc{'bvstrcwmax'}(:,3,3);
rho0 = ncread(url,'rho0');
bstru = squeeze(ncread(url,'bustrcwmax',[i j 1],[1 1 Inf]));
bstrv = squeeze(ncread(url,'bvstrcwmax',[i j 1],[1 1 Inf]));
bstr=abs(complex(bstru,bstrv));
ustr=sqrt(bstr/rho0);
%%
vonKar=0.41;
% gls=nc{'gls'}(:,:,3,3);
% tke=nc{'tke'}(:,:,3,3);
gls = squeeze(ncread(url,'gls',[i j 1 1],[1 1 Inf, Inf]));
tke = squeeze(ncread(url,'tke',[i j 1 1],[1 1 Inf, Inf]));
% gls_av(:,1)=gls(:,1);
% tke_av(:,1)=tke(:,1);
%%
for k=2:nz
gls_av(k,:)=.5*gls(k,:)+.5*gls(k+1,:);
tke_av(k,:)=.5*tke(k,:)+.5*tke(k+1,:);
end
%%
gls_p = -1 ;% ! k-omega
gls_m = 0.5;
gls_n = -1.0;
gls_cmu0 = 0.5477;
exp1 = 3.0+gls_p/gls_n;
exp2 = 1.5+gls_m/gls_n;
exp3 = -1.0/gls_n;
diss0 = (gls_cmu0^exp1).*(tke_av.^exp2).*(gls_av.^exp3);
diss=diss0;
%% diss(:,1)=(bstr.^1.5)/(vonKar*0.9);
if cas<72
diss(1,:)=(ustr.^3)'/(vonKar*0.9);
else
effecz=ustr.*ustar_av.Tr(1:length(ustr))/(2*pi);
diss(1,:)=(ustr.^3)'./(vonKar.*effecz');
end
%%
nu0=1.5e-6;
G0=sqrt(diss0/nu0);
G=sqrt(diss/nu0);
if(0)
for it=1:nt
Gf(:,it)=max([G0(:,it);G(:,it)]);
end
end
%%
% TODO - which G is this...which do we want?
zz = h+z_r;
figure(1); clf
scatter(G(:).^(-.5),1e6*fdiamt(:),14,zz(:),'filled')
set(gca,'fontsize',14)
xlim([0 5])
colorbar('east')
ylabel('Floc Diameter (\mum)')
xlabel('{\itG}^{-1/2}')
%%
s2d = 1. /(3600.*24.);
figure(2);clf
subplot 321
pcolorjw(s2d*tz,z_r,G0)
caxis([0,4])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('G from GLS')
subplot 322
pcolorjw(s2d*tz,z_r,G0)
set(gca,'YLim',[-11.9,-10])
caxis([0,40])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('G from GLS zoom ')
subplot 323
pcolorjw(s2d*tz,z_r,G)
caxis([0,4])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('u^* G ')
subplot 324
pcolorjw(s2d*tz,z_r,G)
set(gca,'YLim',[-11.9,-10])
caxis([0,40])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('u^* G zoom ')
subplot 325
pcolorjw(s2d*tz,z_r,Gf)
caxis([0,4])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('enhanced G ')
subplot 326
pcolorjw(s2d*tz,z_r,Gf)
set(gca,'YLim',[-11.9,-10])
caxis([0,40])
datetick('x','mm/dd','keepticks','keeplimits')
colorbar
title('enhanced G zoom ')
eval(['print -dpng -painters G_',cas,'.png'])