forked from mnemoniko/roms_matlab_sm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sta_mld_sm.m
executable file
·109 lines (94 loc) · 3.34 KB
/
sta_mld_sm.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
function [ mld ] = sta_mld_sm( run )
%STA_MLD_SM This function calculates mixed layer depth for ROMS Ross Sea
%station files based on a threshold method
% INPUT: run - string indicating which run number to analyze
% OUTPUT: mld - mixed layer depth for that run at all stations.
file = romsinitialize_sm('sta',run);
%Basic data
Ipos = nc_varget(file,'Ipos');
Jpos = nc_varget(file,'Jpos');
lat = nc_varget(file,'lat_rho');
lon = nc_varget(file,'lon_rho');
staTime = nc_varget(file,'ocean_time');
staTime = staTime./3600; staTime = staTime./24; %days
staTime = staTime -2190; %set 0 to sept 15 2010
staTime = datenum(2010,9,15+staTime,0,0,0);
%To calc MLDs:
temp = nc_varget(file,'temp');
salt = nc_varget(file,'salt');
zeta = nc_varget(file,'zeta');
h = nc_varget(file, 'h'); %meters
%Remove 0 points:
good = find(Ipos~=0);
Ipos = Ipos(good); Jpos = Jpos(good); temp = temp(:,good,:);
salt = salt(:,good,:); zeta = zeta(:,good);
lat = lat(good); lon = lon(good); h = h(good);
%Set depths of layers
Vtransform = nc_varget(file,'Vtransform');
Vstretching = nc_varget(file,'Vstretching');
theta_s = nc_varget(file,'theta_s');
theta_b = nc_varget(file,'theta_b');
hc = nc_varget(file,'hc'); %critical depth
s_rho = nc_varget(file,'s_rho');
N = length(s_rho);
igrid = 1;
zeta(find(isnan(zeta)))=0; %#ok<*FNDSB>
steps = size(zeta,1); sta = size(zeta,2);
z = zeros(steps,sta,N);
for i=1:steps
for j=1:sta
temp2 = set_depth_sta_sm(Vtransform, Vstretching, theta_s, theta_b, hc, N, igrid, h(j), zeta(i,j),0);
z(i,j,:)=temp2;
end
%display(['Finished depth i = ' num2str(i)]);
end
clear Vtransform Vstretching theta_s theta_b hc N igrid s_rho X Y steps temp2;
display(['Finished setting depth for run ' run]);
rho = zeros(size(z));
for n=1:size(z,1)
for i=1:24
pres = gsw_p_from_z(squeeze(z(n,:,i)),lat);
SA = gsw_SA_from_SP(squeeze(salt(n,:,i)),pres,lon,lat);
CT = gsw_CT_from_pt(SA,squeeze(temp(n,:,i)));
rho(n,:,i) = gsw_rho_CT_exact(SA,CT,10.1325);
end
end
display(['Finished setting rho for run ' run]);
%Set surface layer as reference layer
tempRef = temp(:,:,24);
rhoRef = rho(:,:,24);
%Calculate |deltaT| from surface layer and find first point greater than
%0.2 deg C
%Calculate deltaRho from surface and find first point greater than 0.03
%kg/m^3
deltaT = zeros(size(temp)); deltaRho = deltaT;
for i=1:24
deltaT(:,:,i) = abs(temp(:,:,i)-tempRef);
deltaRho(:,:,i) = rho(:,:,i)-rhoRef;
end
% deltaT(:,:,24)=0; deltaRho(:,:,24)=0;
mldR = zeros(length(staTime),length(Ipos)); mldT = mldR;
for n = 1:length(staTime)
for i=1:length(Ipos)
dT = squeeze(deltaT(n,i,:));
layer = find(dT>.2,1,'last');
if(isempty(layer))
mldT(n,i)=squeeze(z(n,i,1)); %Set to bottom depth
else
mldT(n,i)=squeeze(z(n,i,layer));
end
dR = squeeze(deltaRho(n,i,:));
layer = find(dR>0.03,1,'last');
if(isempty(layer))
mldR(n,i)=squeeze(z(n,i,1));
else
mldR(n,i)=squeeze(z(n,i,layer));
end
end
end
%Compare two depths, and choose shallower of two for fully mixed layer
mld = max(mldT,mldR); %max, because z is negative
display(['Finished calculating mld for run ' run]);
display('Saving data...');
save([run 'MLDsta.mat'],'mld', 'mldR','mldT','run','staTime');
end