/[MITgcm]/MITgcm_contrib/gmaze_pv/compute_MLD.m
ViewVC logotype

Annotation of /MITgcm_contrib/gmaze_pv/compute_MLD.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Fri Oct 6 21:44:53 2006 UTC (17 years, 7 months ago) by gmaze
Branch: MAIN
Changes since 1.1: +11 -1 lines
Ensure MLD > 0

1 gmaze 1.1 %
2     % [] = compute_MLD(SNAPSHOT)
3     %
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     function compute_MLD(snapshot)
32    
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.09998*9.81*Tdpt(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    
159    

  ViewVC Help
Powered by ViewVC 1.1.22