Skip to content

Commit

Permalink
Update for completeness
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhihua-zheng
zhihua-zheng committed Dec 23, 2020
1 parent fb7401f commit ab1489b
Show file tree
Hide file tree
Showing 26 changed files with 823 additions and 352 deletions.
10 changes: 5 additions & 5 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[![DOI](https://zenodo.org/badge/287464756.svg)](https://zenodo.org/badge/latestdoi/287464756)

**This repository contains the MATLAB scripts used in the paper "Evaluating Monin-Obukhov scaling in the unstable oceanic surface layer"**
**This repository contains the MATLAB scripts used in the paper "Evaluating Monin-Obukhov Scaling in the Unstable Oceanic Surface Layer", J. Phys. Oceanogr., doi:10.1175/JPO-D-20-0201.1.**

You can reach me at [[email protected]](mailto:[email protected]) for questions related to the usage of this code.

Expand All @@ -24,7 +24,7 @@ It also uses some external toolboxes/functions, as listed below:
- [tight_subplot](https://www.mathworks.com/matlabcentral/fileexchange/27991-tight_subplot-nh-nw-gap-marg_h-marg_w)

### Data
The data used in this work are mostly publicly available (see the "Data availability statement" in our paper). For convience, a copy of the original data can be downloaded either from [this Google Drive link](https://drive.google.com/file/d/13UYYOT9AXFufjMw6_wr4-hoNv7M3tT7v/view?usp=sharing), or from the most recent [release](https://github.com/zhihua-zheng/EvaMO_UOSL_code/releases/tag/v1.0) of this repository.
The data used in this work are mostly publicly available (see the "Data availability statement" in our paper). For convience, a copy of the original data can be downloaded either from [this Google Drive link](https://drive.google.com/file/d/13UYYOT9AXFufjMw6_wr4-hoNv7M3tT7v/view?usp=sharing), or from the most recent [release](https://github.com/zhihua-zheng/EvaMO_UOSL_code/releases/tag/v1.1) of this repository.

### Work flow
1. Download this repository to your local computer. Alternatively, you can clone it using the command line in Terminal,
Expand All @@ -38,11 +38,11 @@ Add this repository to your MATLAB search path
addpath(genpath('<path to this repository>'))
```

2. Download the data file `EvaMO_UOSL_data.zip` as described in the previous [Data](#-data) section. Unzip it into the `EvaMO_UOSL_code` directory.
2. Download the data file `EvaMO_UOSL_data.zip` following the [Data](#-data) section. Unzip it into the top-level directory of the repository `EvaMO_UOSL_code`.

3. Run the main analysis script `evaMO_main.m`.

### Attribution
This code is freely available for reuse as described in the MIT License. However, if you use this code in an academic publication, a propriate citation to the source would be appreciated:
This set of scripts is free for reuse as described in the MIT License. However, if you use this code in an academic publication, a propriate citation to the source would be appreciated:

* Zhihua Zheng. (2020, August 17). Analysis scripts for "Evaluating Monin-Obukhov scaling in the unstable oceanic surface layer" (Version v1.0). Zenodo. doi:10.5281/zenodo.3988503
* Zhihua Zheng. (2020, December 22). Analysis scripts for "Evaluating Monin-Obukhov Scaling in the Unstable Oceanic Surface Layer". Zenodo. doi:10.5281/zenodo.3988503
26 changes: 14 additions & 12 deletions evaMO_main.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,28 +50,30 @@
plot_profcom(rpPA,1,dTax(2));
plot_profcom(rpSP,2,dTax(4));

%% Figure 5: boxplot for observed phi_h in zeta space
%% Figure 5: boxplot for the observed phi_h in zeta space

phi_obs = {pzPA.phi_qs(:) pzSP.phi_qs(:)};
zet_obs = {pzPA.zeta_qs(:) pzSP.zeta_qs(:)};
PZ_scat_obs(phi_obs,zet_obs);

%% Figure 6: predictions from the 'SE' surface wave breaking model
%% Figure 6: predictions from the "SE" surface wave breaking model

Poff = 0.08;
WBf = figure('position',[0 0 935 530]);
[WBax,WBpos] = tight_subplot(1,2,[.02 .05],[.09 .1],[.05 .09]);
Poff = 0.14;
WBf = figure('position',[0 0 1015 500]);
[WBax,WBpos] = tight_subplot(1,3,[.03 .04],[.09 .1],[.04 .04]);
WBax(1).Position = [WBpos{1}(1:2) WBpos{1}(3)+Poff WBpos{1}(4)];
WBax(2).Position = [WBpos{2}(1)+Poff WBpos{2}(2) WBpos{2}(3)-Poff WBpos{2}(4)];
WBax(2).Position = [WBpos{2}(1)+Poff WBpos{2}(2) WBpos{2}(3)-Poff/2 WBpos{2}(4)];
WBax(3).Position = [WBpos{3}(1)+Poff/2 WBpos{3}(2) WBpos{3}(3)-Poff/2 WBpos{3}(4)];

r2o = -0.2; % z_0/L
plot_SFWB_phi2d(WBax(1),r2o,pzPA,pzSP)
q_in_SFWB(r2o,[],WBf,WBax(2))
q_in_SFWB(WBf,WBax(2),r2o,[])
demo_sfwbProf(WBax(3),r2o)

%% Figure 7: predictions from the SE Langmuir turbulence model

figure('position',[0 0 1015 500]);
[LTax,~] = tight_subplot(1,2,[.03 .03],[.1 .1],[.06 .08]);
[LTax,~] = tight_subplot(1,2,[.03 .03],[.1 .1],[.06 .06]);
plot_H15se_phi2d(LTax,'H2015',pzPA,pzSP);

%% Figure 8: comparison of phi_h from Kansas, observations and models
Expand All @@ -82,14 +84,14 @@
PZ_scat(pzPA,1,PZf,PZax(1),1);
PZ_scat(pzSP,2,PZf,PZax(2),1);

%% Figure 9: phi_h & l^* from the SE Langmuir turbulence model including length scale Eq.
%% Figure 9: phi_h & l^* from the SE Langmuir model with dynamic length scale

% Find optimal E6 value for SMCLT with varying l^*

% dT_obs_SMCLT;
% figure;
% plot(E6_list,slp,'-o'); grid on
% optimal E6 is 2.4, with slp = 1
% optimal E6 is about 2.5, with slp ~ 1
% save([spursi_dir,'spursi_slpE6.mat'],'E6_list','slp')

Lstarf = figure('position',[0 0 820 800]);
Expand All @@ -102,15 +104,15 @@
%% Forcing distribution at two sites

xString = 'H_s [m]';
plot_pdf(skfPA.SKF.SD.Hs(PAqs),skfSP.SKF.SD.Hs(SPqs),22,xString)
plot_pdf(skfPA.SKF.Hs(PAqs),skfSP.SKF.Hs(SPqs),22,xString)
xlim([0 10])

xString = 'u_* [m s^{-1}]';
plot_pdf(skfPA.SKF.Ustar(PAqs),skfSP.SKF.Ustar(SPqs),22,xString)
xlim([0 0.04])

xString = 'B_0 [m^2 s^{-3}]';
plot_pdf(skfPA.SKF.Bo(PAqs),skfSP.SKF.Bo(SPqs),42,xString)
plot_pdf(skfPA.SKF.B0(PAqs),skfSP.SKF.B0(SPqs),42,xString)
xlim([-6e-7 4e-7])

xString = '\eta^x';
Expand Down
96 changes: 96 additions & 0 deletions library/Fin_from_dws.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
function Fin = Fin_from_dws(WQ,fctr,bw,ws,wdir,fop)
%
% Fin_from_dws
%==========================================================================
%
% USAGE:
% Fin = Fin_from_dws(WQ,fctr,bw,ws,wdir,fop)
%
% DESCRIPTION:
% Compute the rate of wind energy input accroding to Plant 1982, or
% Donelan and Pierson 1987.
%
% INPUT:
%
% WQ - timetable for spectral wave data
% fctr - bin center frequency
% bw - bandwidth for frequency bins
% ws - [1,t] wind stress, or [f,t] wind speed at height pi/k
% wdir - [1,t] direction of wind, clockwise from N [degree, from]
% fop - formula option:
% 'P82': Plant 1982
% 'DP87': Donelan and Pierson 1987
%
% OUTPUT:
%
% Fin - the rate of energy input from wind
%
% AUTHOR:
% Sep 6 2020, Zhihua Zheng [ [email protected] ]
%==========================================================================

%% Loading

datm = WQ.datm;
a0 = WQ.wave_spec'/pi;
a1 = WQ.wave_a1';
b1 = WQ.wave_b1';
a2 = WQ.wave_a2';
b2 = WQ.wave_b2';

%% Constants

rhow = 1025;
rhoa = 1.225;
g = 9.81;
fc = fctr(end) + bw(end)/2; % right edge cutoff frequency [Hz]
ntm = length(datm);

%% Reshaping

df = repmat(bw, 1,ntm);
f = repmat(fctr,1,ntm);

%% Resolved spectrum

switch fop
case 'P82'
multi = 0.04*8*pi^3/g/rhoa.*ws;
mu = pi*(cosd(wdir).*a1 + sind(wdir).*b1);
Fin_re = multi .* sum(f.^3 .* mu .* df);

case 'DP87'
c = g/2/pi./f;
multi = 0.194*8*pi^3*rhoa/rhow/g;
mu = pi*(ws.^2.*(a0 + cosd(2*wdir).*a2 + sind(2*wdir).*b2)/2 + ...
a0.*c.^2 - ...
2*ws.*c.*(cosd(wdir).*a1 + sind(wdir).*b1));
Fin_re = multi * sum(f.^3 .* mu .* df);
end

%% Append analytical high frequency tail (Breivik et al. 2014)

a1_tail = a1(end,:);
b1_tail = b1(end,:);

tailC = multi * fc^4;

switch fop
case 'P82'
muC = pi*(cosd(wdir).*a1_tail + sind(wdir).*b1_tail);
Fin_tail = tailC .* muC;

case 'DP87'
muC = pi*(ws(end,:).^2.*(a0(end,:) + cosd(2*wdir).*a2(end,:) + ...
sind(2*wdir).*b2(end,:))/2 + ...
a0(end,:)*c(end)^2 - ...
2*ws(end,:)*c(end).*(cosd(wdir).*a1_tail + ...
sind(wdir).*b1_tail));
Fin_tail = tailC * muC;
end

%% Add up

Fin = Fin_re + Fin_tail;

end
20 changes: 10 additions & 10 deletions library/MOSTpar_from_flux.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@
waterType = 3;
end

w_b_0 = SKF.w_b_0';
w_s_0 = SKF.w_s_0';
w_theta_0 = SKF.w_theta_0';
alpha = SKF.alpha';
NSW = SKF.nsw;
w_b_0 = SKF.w_b_0';
w_s_0 = SKF.w_s_0';
w_t_0 = SKF.w_t_0';
alpha = SKF.alpha';
NSW = SKF.nsw;

%% Constants

Expand Down Expand Up @@ -81,8 +81,8 @@
end

% depth dependent flux due to the penetrative shortwave radiation
w_theta_r = -(NSW' - Iz)/cp/rho0;
w_b_r = g*(alpha .* w_theta_r);
w_t_r = -(NSW' - Iz)/cp/rho0;
w_b_r = g*(alpha .* w_t_r);

%% MOST parameters

Expand All @@ -92,7 +92,7 @@
% surface layer scales
FS.Ustar = SKF.Ustar';
FS.Sstar = -w_s_0 ./ FS.Ustar / kappa;
FS.Tstar = -(w_theta_r + w_theta_0) ./ FS.Ustar / kappa;
FS.Tstar = -(w_t_r + w_t_0) ./ FS.Ustar / kappa;
FS.Bstar = FS.Bf ./ FS.Ustar / kappa;

% Monin-Obukhov length [m]
Expand All @@ -106,7 +106,7 @@
end

% Append non-turbuelnt fluxes
FS.w_theta_r = w_theta_r;
FS.w_b_r = w_b_r;
FS.w_t_r = w_t_r;
FS.w_b_r = w_b_r;

end
31 changes: 21 additions & 10 deletions library/PZ_scat.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,12 @@ function PZ_scat(pzData,iD,fh,axm,fpart)

zet_range(1,:) = [2e-3 5];
zet_range(2,:) = [3e-2 5];
% zet_range(2,:) = [2e-3 5];
zet_range(3,:) = [3e-3 5];

bin_lim(1,:) = [-3.0 0.7];
bin_lim(2,:) = [-1.6 0.7];
% bin_lim(2,:) = [-3.0 0.7];
bin_lim(3,:) = [-3.0 0.7];

phi_range(1,:) = [6e-2 1.5];
Expand Down Expand Up @@ -36,19 +38,23 @@ function PZ_scat(pzData,iD,fh,axm,fpart)
etaY = pzData.etaY_qs;
fzS = pzData.fzS_qs;
xi = pzData.xi_qs;
alB = pzData.alB_qs;
alB = repmat(alB,size(phi,1),1);

%% Preprocessing

xi(xi<1) = NaN;
xi(xi<1) = NaN;
alB(alB<0) = NaN;
Igood = ~isnan(phi) & ~isnan(zeta) & ~isnan(etaX) &...
~isnan(etaY) & ~isnan(fzS) & ~isnan(xi);
~isnan(etaY) & ~isnan(fzS) & ~isnan(xi);...& ~isnan(alB);

x_data = zeta(Igood);
y_data = phi(Igood);
etX_data = etaX(Igood);
etY_data = etaY(Igood);
fzS_data = fzS(Igood);
xi_data = xi(Igood);
alB_data = alB(Igood);

IposX = x_data > 0;
InegX = x_data <= 0;
Expand All @@ -59,20 +65,22 @@ function PZ_scat(pzData,iD,fh,axm,fpart)
pey_data = etY_data(IposX);
pf_data = fzS_data(IposX);
pxi_data = xi_data(IposX);
pal_data = alB_data(IposX);

ny_data = y_data(InegX);
nx_data = -x_data(InegX);
ny_data = y_data(InegX);
nx_data = -x_data(InegX);
nex_data = etX_data(InegX);
ney_data = etY_data(InegX);
nf_data = fzS_data(InegX);
nf_data = fzS_data(InegX);
nxi_data = xi_data(InegX);
nal_data = alB_data(InegX);

%% Figure

switch fpart
case 1 % Unstable side

Ifrcc = nx_data <= 3; % "forced convection"
Ifrcc = nx_data <= 3; % forced convection
bin_xi = logspace(bin_lim(iD,1),bin_lim(iD,2),24)';
Ngood = numel(nx_data(Ifrcc));
NminC = Ngood*0.005;
Expand All @@ -84,7 +92,10 @@ function PZ_scat(pzData,iD,fh,axm,fpart)
Nphi_H15 = get_H15se_phi_nol(-nx_data,nex_data,ney_data,nf_data,2);
Nphi_H15B = get_H15se_phi_nol(-nx_data,nex_data,ney_data,0,2);
Nphi_H15T = get_H15se_phi_nol(-nx_data,nex_data,ney_data,1,2);
Nphi_WBv1 = get_SFWB_phi_apr(-nx_data,nxi_data,1);

constAlB = repmat(100,size(nal_data));
Nphi_WBv1 = get_SFWB_phi_apr(-nx_data,nxi_data,constAlB,1);
% Nphi_WBv1 = get_SFWB_phi_apr(-nx_data,nxi_data,nal_data,1);

Nref71 = plot(axm,-Nzet,Nphi71,'color',[.5 .5 .5],...
'linewidth',3);
Expand Down Expand Up @@ -198,9 +209,9 @@ function PZ_scat(pzData,iD,fh,axm,fpart)

Pzet = logspace(-5,1.5,100);
[Pphi71,~] = get_emp_phi(Pzet,'Kansas');
Pphi_H15 = get_H15se_phi_nol(px_data,pex_data,pey_data,pf_data,2);
Pphi_H15B = get_H15se_phi_nol(px_data,pex_data,pey_data,0,2);
Pphi_H15T = get_H15se_phi_nol(px_data,pex_data,pey_data,1,2);
Pphi_H15 = get_H15se_phi_nol(px_data,pex_data,pey_data,pf_data,2);
Pphi_H15B = get_H15se_phi_nol(px_data,pex_data,pey_data,0,2);
Pphi_H15T = get_H15se_phi_nol(px_data,pex_data,pey_data,1,2);

Pgood = numel(px_data);
PminC = max(Pgood*0.001,5);
Expand Down
Loading

0 comments on commit ab1489b

Please sign in to comment.