Change to relative directory path

.gitignore

@@ -0,0 +1,2 @@
+# Ignore data directory
+/data

README.md

@@ -5,30 +5,46 @@
 You can reach me at [[email protected]](mailto:[email protected]) for questions related to the usage of this code. 
 
 ### Description
-This set of scripts fufills three stages in generating the results presented in the paper:
-- Load and process raw data
+This set of scripts fufills the analysis presented in the paper in three steps:
+- Load and process original data
 - Calculate key parameters
 - Produce figures
 
 ### Prerequisites
 This code is tested in MATLAB R2019b.
 
 It also uses some external toolboxes/functions, as listed below:
-- [Gibbs-SeaWater (GSW) Oceanographic Toolbox'](http://www.teos-10.org/software.htm)
+- [Gibbs-SeaWater (GSW) Oceanographic Toolbox](http://www.teos-10.org/software.htm)
 - [Air-sea toolbox](https://github.com/sea-mat/air-sea)
 - [cbrewer](https://www.mathworks.com/matlabcentral/fileexchange/34087-cbrewer-colorbrewer-schemes-for-matlab)
-- [RGB](https://www.mathworks.com/matlabcentral/fileexchange/24497-rgb-triple-of-color-name-version-2) triple of color name
+- [Custom Colormap](https://www.mathworks.com/matlabcentral/fileexchange/69470-custom-colormap)
+- [RGB](https://www.mathworks.com/matlabcentral/fileexchange/24497-rgb-triple-of-color-name-version-2)
+- [Marke​rTransparency](https://www.mathworks.com/matlabcentral/fileexchange/65194-peterrochford-markertransparency)
+- [suplabel](https://www.mathworks.com/matlabcentral/fileexchange/7772-suplabel)
+- [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 your 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]() of this repository.
 
 ### Work flow
+1. Download this repository to your local computer. Alternatively, you can clone it using the command line in Terminal,
 
-### Reference
+```bash
+git clone https://github.com/zhihua-zheng/EvaMO_UOSL_code.git
+```
+Add this repository to your MATLAB search path
+
+```matlab
+addpath(genpath('<path to this repository>'))
+```
+
+2. Download the data file `EvaMO_UOSL_data.zip` as described in the [Data](#-data) section. Unzip it into the `EvaMO_UOSL_code` directory and rename the `EvaMO_UOSL_data` direcctory to `data`.
+
+3. 
+
+4.
 
-* McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.
 
-* Charles, R. 2015: cbrewer : colorbrewer schemes for Matlab, MATLAB Central File Exchange. Retrieved May 26, 2019.
 
 Functions for solving super-equilibirum model equations are in ./library
 ========================================================================

evaMO_main.m

@@ -1,9 +1,25 @@
+%% evaMO_main
+
+%% Read original data
+
+read_ocsp_TS(0);
+read_ocsp_flux(0);
+read_ocsp_wave(0);
+
+read_spursi_TS;
+read_spursi_flux;
+read_spursi_wave;
+
+%% Process original data
+
+prep_ocsp;
+prep_spursi;
 
 %% Load processed datasets
 
-root_dir   = '~/GDrive/UW/Research/Data/';
-ocsp_dir   = [root_dir,'Papa/'];
-spursi_dir = [root_dir,'SPURSI/'];
+data_dir   = './data/';
+ocsp_dir   = [data_dir,'Papa/'];
+spursi_dir = [data_dir,'SPURSI/'];
 
 pzPA = load([ocsp_dir,  'ocsp_pzData.mat']);
 pzSP = load([spursi_dir,'spursi_pzData.mat']);
@@ -74,11 +90,12 @@
 % figure;
 % plot(E6_list,slp,'-o'); grid on
 % optimal E6 is 2.4, with slp = 1
-save([spursi_dir,'spursi_slpE6.mat'],'E6_list','slp')
+% save([spursi_dir,'spursi_slpE6.mat'],'E6_list','slp')
 
 Lstarf = figure('position',[0 0 820 800]);
 [lsrAx,~] = tight_subplot(2,2,[.07 .02],[.07 .03],[.08 .03]);
 
+oE6 = 2.5;
 phi_lstar_E6(pzPA,1,oE6,lsrAx,1);
 phi_lstar_E6(pzSP,2,oE6,lsrAx,1);
 
@@ -88,11 +105,11 @@
 plot_pdf(skfPA.SKF.SD.Hs(PAqs),skfSP.SKF.SD.Hs(SPqs),22,xString)
 xlim([0 10])
 
-xString = 'u_* [m/s]';
+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]';
+xString = 'B_0 [m^2 s^{-3}]';
 plot_pdf(skfPA.SKF.Bo(PAqs),skfSP.SKF.Bo(SPqs),42,xString)
 xlim([-6e-7 4e-7])
 

library/MOSTpar_from_flux.m

@@ -12,16 +12,17 @@
 %
 % INPUT:
 %
-%  idatm - 1-D column vector of MATLAB datetime for the inquired period
-%  zz - vertical coordinates for the choosen 2 levels [-, m]
-%  casename - string to indicate the parameters for which cases is inquired
+%  idatm    - 1-D column vector, MATLAB datetime for the period inquired
+%  zz       - Vertical coordinates for the choosen 2 levels [-, m]
+%  casename - Name of the case for which parameters are inquired
 %         
 % OUTPUT:
 %
-%  MOL - Monin-Obukhov length for evaluated depth [m]
+%  MOL  - Monin-Obukhov length for evaluated depth [m]
 %  zeta - Monin-Obukhov stability parameter
-%  FS - struct contains turbulence scale for velocity (Ustar), temperature 
-%       (Tstar), salinity (Sstar), buoyancy (Bstar), buoyancy forcing Bf
+%  FS   - Struct contains turbulence scale for velocity (Ustar), 
+%         temperature (Tstar), salinity (Sstar), buoyancy (Bstar),
+%          buoyancy forcing Bf
 %
 % AUTHOR:
 %  June 29 2019, Zhihua Zheng                             [ [email protected] ]
@@ -48,22 +49,10 @@
 
 ntm = length(w_b_0);
 
-MOconsts_name = '~/GDrive/UW/Research/Data/Misce/MO_consts.mat';
-load(MOconsts_name,'g','rho0','cp','kappa');
-
-%% T-S dependent coefficients
-
-% ignore the depth-dependence as T-S variation with depth is small
-
-% use sea surface SA and CT before calibration
-% ssSA = gsw_SA_from_SP(SKF.sss,ssp,lon,lat); % [g/kg]
-% ssCT = gsw_CT_from_t(ssSA,SKF.sst,ssp);
-% 
-% % isobaric heat capacity of seawater [J/kg/C]
-% % cp    = gsw_cp_t_exact(ssSA,SKF.sst,ssp);
-% 
-% alpha = gsw_alpha(ssSA,ssCT,ssp);
-% % beta  = gsw_beta(ssSA,ssCT,ssp);
+g     = 9.81;
+rho0  = 1025;
+cp    = 3985;
+kappa = 0.4;
 
 %% Fluxes due to shortwave radiation
 
@@ -81,13 +70,13 @@
 
 switch band_SR
 
-    case 0     % test: shortwave radiation are absorbed at the surface
+    case 0 % test: shortwave radiation are absorbed at the surface
     Iz = zeros(ndum,ntm);
 
-    case 2     % 2-band exponential decay of shortwave radiation
+    case 2 % 2-band exponential decay of shortwave radiation
     Iz = get_SRz(NSW,zdum,band_SR,waterType);
 
-    case 9     % 9-band exponential decay of shortwave radiation
+    case 9 % 9-band exponential decay of shortwave radiation
     Iz = get_SRz(NSW,zdum,band_SR,waterType);
 end
 

library/get_MOST_prof.m

@@ -33,22 +33,6 @@
 %==========================================================================
 
 [~,ntm] = size(zeta);
-
-% Use fine resolution for integration between each level of z
-% This integration assumes depth independent Tstar, MOL
-% MOL = abs(z) ./ zeta;
-% cdPTmost = zeros(nz,ntm); % accumulative dPT
-% for i = 2:nz
-%     
-%     iz = linspace(z(1),z(i))';
-%     izeta = abs(iz) ./ repmat(MOL(1,:),100,1);
-%     iTstar = repmat(Tstar(1,:),100,1);
-%     phi_dum = get_emp_phi(izeta,fop);
-%     
-%     intgrd = iTstar ./ abs(iz) .* phi_dum;
-%     cdPTmost(i,:) = trapz(iz,intgrd);
-% end
-
 phi_dum = get_emp_phi(zeta,fop);
 intgrd  = Tstar ./ abs(zdum) .* phi_dum;
 

library/get_Rib.m

@@ -26,9 +26,7 @@
 
 %% Presetting
 
-MOconsts_name = '~/GDrive/UW/Research/Data/Misce/MO_consts.mat';
-load(MOconsts_name,'kappa');
-
+kappa = 0.4;
 nz  = length(z);
 ntm = length(Ustar);
 
@@ -48,5 +46,4 @@
 
 Ribprof = N2 ./ S2;
 
-end
-
+end

library/get_SMCse_phi_nol.m

@@ -84,7 +84,7 @@
 phiS = nan(size(zeta));
 
 fun = @SMCse_phiEq_nol;
-options = optimovtions('fsolve',...'Algorithm','Levenberg-Marquardt',...
+options = optimoptions('fsolve',...'Algorithm','Levenberg-Marquardt',...
                        'Display','off');
 
 for i = 1:nzet

library/get_fit_gdT.m

@@ -1,4 +1,4 @@
-function [gradT,PTfit] = get_fit_gdT(depth,WLD,PTprof,BLD,nFit,ishow)
+function [gradT,PTfit] = get_fit_gdT(depth,iD,PTprof,BLD,nFit,ishow)
 
 %% Pre-setting
 
@@ -25,11 +25,10 @@
     % number of measurements within surface layer
     iSLDj   = sum(depth <= SLD(j));
     nPTinSL = sum(~isnan(PTprof(1:mlim(2),j)));
-    iWLDj   = sum(depth <= WLD(j))+1;
 
     %% loop for regressions with varying points
 
-    if nPTinSL >= nFit && ~isnan(WLD(j))
+    if nPTinSL >= nFit
 
         p   = cell(1,nr);
         gof = cell(1,nr);
@@ -41,7 +40,7 @@
 
         for ir = 1:nr
 
-        Ifit = ~isnan(PTprof(1:m(ir),j));% & depth(1:m(ir)) > WLD(j);
+        Ifit = ~isnan(PTprof(1:m(ir),j));
 
         if sum(Ifit) >= 3
         [p{ir},gof{ir}] = fit(lnz(Ifit),PTprof(Ifit,j),'poly2');
@@ -87,17 +86,21 @@
 
             % close to SL base, gradient is small, susceptible to error
             nGrad   = min([M-1 iSLDj]);
-            Ieva    = depth(1:nGrad) > depth(1);% WLD(j);
+            Ieva    = depth(1:nGrad) > depth(1);
             d_eva   = depth(Ieva);
             lnz_eva = log(d_eva);
             gradT(Ieva,j) = -(2*b(1)*lnz_eva + b(2)) ./ d_eva;
 
             % fitted profile in surface layer
-            PTfit(iWLDj:iSLDj,j) = P(lnz(iWLDj:iSLDj));
+            PTfit(1:iSLDj,j) = P(lnz(1:iSLDj));
 
             % visualize polynomial regression
             if ishow
                 show_polyfit_lnz;
+
+            % Figure 3 is from the SPURS-I dataset with j = 3427
+            elseif iD == 2 && j == 3427
+                show_polyfit_lnz;
             end
         end
 
@@ -152,8 +155,6 @@ function show_polyfit_lnz()
     ylh.Units = 'Normalized';
     set(xlh,'Position',xlh.Position + [0 -off_r 0]);
     set(ylh,'Position',ylh.Position + [off_r 0 0]);
-
-% demo_prof_fit_SP.png is from j = 3427
 end
 
 end

library/get_qs_time.m

@@ -43,6 +43,8 @@
 
 %% Show quasi-steady state periods
 
+if iD == 2
+
 figure('position',[0 0 900 600]);
 [hax,~] = tight_subplot(2,1,[.03 .02],[.1 .05],[.1 .1]);
 
@@ -101,3 +103,4 @@
        'location','north')
 
 set(hax,'TickDir','out')
+end

library/get_skf_prof.m

@@ -2,13 +2,13 @@
 
 %% Parsing inputs
 
-root_dir = '~/GDrive/UW/Research/Data/';
+data_dir = './data/';
 cLave = num2str(Lave);
 
 switch casename
 
     case 'Papa'
-        data_dir = [root_dir,'Papa/'];
+        data_dir = [data_dir,'Papa/'];
         PROFname = fullfile(data_dir,['ocsp_prof_',cLave,'hrPMEL.mat']);
         FLUXname = fullfile(data_dir,['ocsp_flux_',cLave,'hrPMEL.mat']);
         WAVEname = fullfile(data_dir,['ocsp_wave_',cLave,'hr.mat']);
@@ -17,7 +17,7 @@
         ssp      = 1;
 
     case 'SPURSI'
-        data_dir = [root_dir,'SPURSI/'];
+        data_dir = [data_dir,'SPURSI/'];
         PROFname = fullfile(data_dir,['spursi_prof_',cLave,'hrBox.mat']);        
         FLUXname = fullfile(data_dir,['spursi_flux_',cLave,'hrUOP.mat']);
         WAVEname = fullfile(data_dir,['spursi_wave_',cLave,'hr.mat']);
@@ -91,8 +91,10 @@
 
 %% Constants
 
-MOconsts_name = [root_dir,'Misce/MO_consts.mat'];
-load(MOconsts_name,'g','rho0','cp','kappa');
+g     = 9.81;
+rho0  = 1025;
+cp    = 3985;
+kappa = 0.4;
 
 %% Seawater expansion/contraction coefficients
 

library/ocsp_rp.m

@@ -1,12 +1,12 @@
-function [rPTprof,rPTmost,zml,Iprof] = ocsp_rp(PTprof,depth_t,...
+function [rPTprof,rPTmost,zbl,Iprof] = ocsp_rp(PTprof,depth_t,...
                                        BLD,dataName,SKF,Islc)
 
 SLD = BLD/5;
 
 % find the mean temperature in surface layer
 mPTprof = get_mT_inLayer(PTprof,depth_t,SLD,@nanmean);
 rPTprof = PTprof - mPTprof; % relative to mPTprof
-zml     = -depth_t ./ BLD';
+zbl     = -depth_t ./ BLD';
 
 % temperature profile from MOST integration
 zdum = -[(1:.1:20) (21:100) (105:5:300)]';
@@ -17,7 +17,7 @@
 
 %% Indices for unstable side
 
-IspanSL = max(zml) > -0.02 & min(zml) < -0.2;
+IspanSL = max(zbl) > -0.02 & min(zbl) < -0.2;
 Ineg    = sum(MOL>0) == 0;
 Iprof   = Islc & Ineg & IspanSL;
 

library/pinBin.m

@@ -96,7 +96,7 @@
 
 % critical vlaues for Student's distribution
 S.tupp = tinv(alpU,S.dof);
-% S.tlow = tinv(alpL,S.dof);
+S.tlow = tinv(alpL,S.dof);
 
 S.qerr = S.tupp .* S.qstd ./ sqrt(S.dof); % symmetric confidence limits
 

library/plot_H15se_phi2d.m

@@ -202,7 +202,7 @@ function plot_H15se_phi2d(LTax,mop,pzPA,pzSP)
       [char(981),'_h(',char(950),', ',char(402),'^{s}_z) with ',...
        char(951),'^y = 0, ',char(951),'^x = ',num2str(cEtX)],...
       [char(981),'_h(',char(950),', ',char(402),'^{s}_z)/',...
-       char(981),'_h^{No wave}'],...
+       char(981),'_h^{no wave}'],...
       'OCSP data','SPURS-I data',...
       'fontsize',16,'location','northwest','Interpreter','tex');
 hl4.Position = [hl4.Position(1) hl4.Position(2)-0.11 hl4.Position(3:4)];

library/plot_pdf.m

@@ -32,7 +32,7 @@ function plot_pdf(xdata1,xdata2,numB,xString)
 ax = gca;
 ax.Position = ax.Position + off_ratio*[1 1 -1 -1];
 ylh = ylabel('Normalized PDF','fontsize',16);
-xlh = xlabel(xString,'fontsize',16);
+xlh = xlabel(xString,'fontsize',16,'Interpreter','tex');
 xlh.Units = 'Normalized';
 ylh.Units = 'Normalized';
 set(xlh,'Position',xlh.Position + [0 -off_ratio/6 0]);