MITgcm_contrib/gmaze_pv/compute_MLD.m

%
% [MLD] = compute_MLD(SNAPSHOT)
%
% Here we compute the Mixed Layer Depth as:
% MLD = min depth for which : ST > ST(SSS,SST-0.8,p0)  
%
% where:
%  ST is potential density (kg/m3)
%  SST the Sea Surface Temperature (oC)
%  SSS the Sea Surface Salinity (PSU-35)
%  p0  the Sea Level Pressure (mb)
%  EKL is the Ekman layer depth (m, positive)
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_SIGMATHETA>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
% OUTPUT
% ./netcdf-files/<SNAPSHOT>/<netcdf_MLD>.<netcdf_domain>.<netcdf_suff>
% 
% with netcdf_* as global variables
% netcdf_MLD = 'MLD' by default
%
% Rq: This method leads to a MLD deeper than KPPmld in the middle of the 
% ocean, and shallower along the coast.
%
% 09/20/06
% gmaze@mit.edu

function varargout = compute_MLD(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_SIGMATHETA netcdf_THETA netcdf_SALTanom netcdf_MLD
pv_checkpath


% NETCDF file name:
filST = netcdf_SIGMATHETA;
filT  = netcdf_THETA;
filS  = netcdf_SALTanom;

% Path and extension to find them:
pathname = strcat('netcdf-files',sla);
ext = netcdf_suff;

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
ncST    = netcdf(ferfile,'nowrite');
ST      = ncST{4}(:,:,:);
[STlon STlat STdpt] = coordfromnc(ncST);

ferfile = strcat(pathname,sla,snapshot,sla,filT,'.',netcdf_domain,'.',ext);
ncT    = netcdf(ferfile,'nowrite');
SST      = ncT{4}(1,:,:);
[Tlon Tlat Tdpt] = coordfromnc(ncT);

ferfile = strcat(pathname,sla,snapshot,sla,filS,'.',netcdf_domain,'.',ext);
ncS   = netcdf(ferfile,'nowrite');
SSS     = ncS{4}(1,:,:);
[Slon Slat Sdpt] = coordfromnc(ncS);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE The Mixed Layer Depth:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow, disp('pre-allocate'), end
nx = length(STlon);
ny = length(STlat);
SST08 = SST - 0.8;
SSS   = SSS + 35;
Surfadens08 = densjmd95(SSS,SST08,(0.09998*9.81*Tdpt(1))*ones(ny,nx))-1000;
MLD = zeros(size(ST,2),size(ST,3));

if toshow, disp('get MLD'), end
for iy = 1 : size(ST,2)
  for ix = 1 : size(ST,3)
      mm =  find( squeeze(ST(:,iy,ix)) > Surfadens08(iy,ix) );
      if ~isempty(mm)
        MLD(iy,ix) = STdpt(min(mm));
      end
    %end
  end
end

MLD(isnan(squeeze(ST(1,:,:)))) = NaN;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ensure we have the right sign (positive)
mm = nanmean(nanmean(MLD,1));
if mm <= 0
  MLD = -MLD;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Record
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow, disp('record'), end

% General informations: 
if ~isempty('netcdf_MLD')
  netfil = netcdf_MLD;
else
  netfil = 'MLD';
end
units      = 'm';
ncid       = 'MLD';
longname   = 'Mixed Layer Depth';
uniquename = 'MLD';

% Open output file:
nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');

% Define axis:
nx = length(STlon) ;
ny = length(STlat) ;
nz = 1 ;

nc('X') = nx;
nc('Y') = ny;
nc('Z') = nz;
 
nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = STlon;
 
nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = STlat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = STdpt(1);
 
% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = MLD;

nc=close(nc);
close(ncST);
close(ncS);
close(ncT);


% Output:
output = struct('MLD',MLD,'lat',DTlat,'lon',STlon);
switch nargout
 case 1
  varargout(1) = {output};
end
1	gmaze	1.1	%
2	gmaze	1.3	% [MLD] = compute_MLD(SNAPSHOT)
3	gmaze	1.1	%
4			% Here we compute the Mixed Layer Depth as:
5			% MLD = min depth for which : ST > ST(SSS,SST-0.8,p0)
6			%
7			% where:
8			% ST is potential density (kg/m3)
9			% SST the Sea Surface Temperature (oC)
10			% SSS the Sea Surface Salinity (PSU-35)
11			% p0 the Sea Level Pressure (mb)
12	gmaze	1.2	% EKL is the Ekman layer depth (m, positive)
13	gmaze	1.1	%
14			% Files names are:
15			% INPUT:
16			% ./netcdf-files/<SNAPSHOT>/<netcdf_SIGMATHETA>.<netcdf_domain>.<netcdf_suff>
17			% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
18			% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
19			% OUTPUT
20			% ./netcdf-files/<SNAPSHOT>/<netcdf_MLD>.<netcdf_domain>.<netcdf_suff>
21			%
22			% with netcdf_* as global variables
23			% netcdf_MLD = 'MLD' by default
24			%
25			% Rq: This method leads to a MLD deeper than KPPmld in the middle of the
26			% ocean, and shallower along the coast.
27			%
28			% 09/20/06
29			% gmaze@mit.edu
30
31	gmaze	1.3	function varargout = compute_MLD(snapshot)
32	gmaze	1.1
33			global sla toshow
34			global netcdf_suff netcdf_domain
35			global netcdf_SIGMATHETA netcdf_THETA netcdf_SALTanom netcdf_MLD
36			pv_checkpath
37
38
39			% NETCDF file name:
40			filST = netcdf_SIGMATHETA;
41			filT = netcdf_THETA;
42			filS = netcdf_SALTanom;
43
44			% Path and extension to find them:
45			pathname = strcat('netcdf-files',sla);
46			ext = netcdf_suff;
47
48			% Load files:
49			ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
50			ncST = netcdf(ferfile,'nowrite');
51			ST = ncST{4}(:,:,:);
52			[STlon STlat STdpt] = coordfromnc(ncST);
53
54			ferfile = strcat(pathname,sla,snapshot,sla,filT,'.',netcdf_domain,'.',ext);
55			ncT = netcdf(ferfile,'nowrite');
56			SST = ncT{4}(1,:,:);
57			[Tlon Tlat Tdpt] = coordfromnc(ncT);
58
59			ferfile = strcat(pathname,sla,snapshot,sla,filS,'.',netcdf_domain,'.',ext);
60			ncS = netcdf(ferfile,'nowrite');
61			SSS = ncS{4}(1,:,:);
62			[Slon Slat Sdpt] = coordfromnc(ncS);
63
64
65			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66			% COMPUTE The Mixed Layer Depth:
67			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68			if toshow, disp('pre-allocate'), end
69			nx = length(STlon);
70			ny = length(STlat);
71			SST08 = SST - 0.8;
72			SSS = SSS + 35;
73			Surfadens08 = densjmd95(SSS,SST08,(0.099989.81Tdpt(1))*ones(ny,nx))-1000;
74			MLD = zeros(size(ST,2),size(ST,3));
75
76			if toshow, disp('get MLD'), end
77			for iy = 1 : size(ST,2)
78			for ix = 1 : size(ST,3)
79			mm = find( squeeze(ST(:,iy,ix)) > Surfadens08(iy,ix) );
80			if ~isempty(mm)
81			MLD(iy,ix) = STdpt(min(mm));
82			end
83			%end
84			end
85			end
86
87			MLD(isnan(squeeze(ST(1,:,:)))) = NaN;
88
89	gmaze	1.2
90			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91			% Ensure we have the right sign (positive)
92			mm = nanmean(nanmean(MLD,1));
93			if mm <= 0
94			MLD = -MLD;
95			end
96
97
98
99	gmaze	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100			% Record
101			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
102			if toshow, disp('record'), end
103
104			% General informations:
105			if ~isempty('netcdf_MLD')
106			netfil = netcdf_MLD;
107			else
108			netfil = 'MLD';
109			end
110			units = 'm';
111			ncid = 'MLD';
112			longname = 'Mixed Layer Depth';
113			uniquename = 'MLD';
114
115			% Open output file:
116			nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
117
118			% Define axis:
119			nx = length(STlon) ;
120			ny = length(STlat) ;
121			nz = 1 ;
122
123			nc('X') = nx;
124			nc('Y') = ny;
125			nc('Z') = nz;
126
127			nc{'X'} = ncfloat('X');
128			nc{'X'}.uniquename = ncchar('X');
129			nc{'X'}.long_name = ncchar('longitude');
130			nc{'X'}.gridtype = nclong(0);
131			nc{'X'}.units = ncchar('degrees_east');
132			nc{'X'}(:) = STlon;
133
134			nc{'Y'} = ncfloat('Y');
135			nc{'Y'}.uniquename = ncchar('Y');
136			nc{'Y'}.long_name = ncchar('latitude');
137			nc{'Y'}.gridtype = nclong(0);
138			nc{'Y'}.units = ncchar('degrees_north');
139			nc{'Y'}(:) = STlat;
140
141			nc{'Z'} = ncfloat('Z');
142			nc{'Z'}.uniquename = ncchar('Z');
143			nc{'Z'}.long_name = ncchar('depth');
144			nc{'Z'}.gridtype = nclong(0);
145			nc{'Z'}.units = ncchar('m');
146			nc{'Z'}(:) = STdpt(1);
147
148			% And main field:
149			nc{ncid} = ncfloat('Z', 'Y', 'X');
150			nc{ncid}.units = ncchar(units);
151			nc{ncid}.missing_value = ncfloat(NaN);
152			nc{ncid}.FillValue_ = ncfloat(NaN);
153			nc{ncid}.longname = ncchar(longname);
154			nc{ncid}.uniquename = ncchar(uniquename);
155			nc{ncid}(:,:,:) = MLD;
156
157			nc=close(nc);
158	gmaze	1.3	close(ncST);
159			close(ncS);
160			close(ncT);
161	gmaze	1.1
162
163	gmaze	1.3	% Output:
164			output = struct('MLD',MLD,'lat',DTlat,'lon',STlon);
165			switch nargout
166			case 1
167			varargout(1) = {output};
168			end