Warning, /verification/tutorial_baroclinic_gyre/analysis/matlab_plots.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit ce0d9af5 on 2021-04-21 19:30:20 UTC
ce0d9af5ea Jeff*0001 % read in additional colormaps for plots
0002 % assumes 'addpath ../../../utils/matlab/' typed from your run directory
0003 % which is also where other matlab functions such as 'rdmds' are located
0004 bluered_colormaps
94151a9b18 Jeff*0005
ce0d9af5ea Jeff*0006
0007 %%%%%%%%%%% load grid data %%%%%%%%%%%
0008
0009 % Load grid variables; names correspond with MITgcm source code.
0010 % (if using standard binary output instead of netcdf, these are dumped
0011 % to individual files, e.g. 'RC.data' etc.)
0012 % Assumes that output in multiple tiles has been concatenated into
0013 % global files (used script utils/python/MITgcmutils/scripts/gluemncbig)
0014 % and moved into the top directory;
0015 % see MITgcm user manual section 4.2.4.1
0016
0017 % Using a spherical polar grid, all X,Y variables in lon,lat coordinates
0018 % vertical grid provided in meters, area in m^2
0019
0020 % 1-D fields
0021 RC = ncread('grid.nc', 'RC'); % vertical grid, cell center locations
0022 drF = ncread('grid.nc', 'drF'); % vertical spacing of grid cells (thickness of cells)
0023 % 2-D fields (x,y)
0024 XC = ncread('grid.nc', 'XC'); % x-location of gridcell centers
0025 YC = ncread('grid.nc', 'YC'); % y-location of gridcell centers
0026 dyG = ncread('grid.nc', 'dyG'); % grid spacing in y-dim (separation between corners)
0027 rA = ncread('grid.nc', 'rA'); % surface area of gridcells
0028 % 3-D fields (x,y,z)
0029 HFacC = ncread('grid.nc', 'HFacC'); % vertical fraction of cell which is ocean
94151a9b18 Jeff*0030
ce0d9af5ea Jeff*0031 % For convenience, load additional dimensional variables from
0032 % netcdf files (these variables are NOT included in binary output)
0033 % These grid variables are also present in netcdf diagnostic data files.
0034 % 1-D fields
0035 X = ncread('grid.nc', 'X'); % 1-D version of XC data
0036 Y = ncread('grid.nc', 'Y'); % 1-D version of YC data
0037 Xp1 = ncread('grid.nc', 'Xp1'); % x-location of lower left corner
0038 Yp1 = ncread('grid.nc', 'Yp1'); % y-location of lower left corner
0039 % Xp1 and Yp1 are effectively 1-D versions of grid variables XG,YG with
0040 % an extra data value at the eastern and northern ends of the domain.
0041 % See MITgcm users manual section 2.11 for additional MITgcm grid info
94151a9b18 Jeff*0042
ce0d9af5ea Jeff*0043 % Number of gridcells in x,y, full domain:
0044 Nx = size(XC, 1);
0045 Ny = size(XC, 2);
0046 % out-of-the-box tutorial setup is configured to be 62x62
94151a9b18 Jeff*0047
ce0d9af5ea Jeff*0048
0049 %%%%%%%%%%% load diagnostics %%%%%%%%%%%
94151a9b18 Jeff*0050
ce0d9af5ea Jeff*0051 % unit for temperature is degrees Celsius, velocity in m/s,
0052 % surface height in m, heat flux in W/m^2
0053 % see run output file available_diagnostics.log
94151a9b18 Jeff*0054
ce0d9af5ea Jeff*0055 % load statistical diagnostic output (set to monthly time-avg output)
0056 % only one output region is defined: global (the default)
0057 TRELAX_ave = ncread('dynStDiag.nc', 'TRELAX_ave'); % (depth, region, time); depth=1 (surf-only), region=1 (global)
0058 THETA_lv_ave = ncread('dynStDiag.nc', 'THETA_lv_ave'); % (depth, region, time); region=1 (global)
0059 THETA_lv_std = ncread('dynStDiag.nc', 'THETA_lv_std'); % (depth, region, time); region=1 (global)
94151a9b18 Jeff*0060
ce0d9af5ea Jeff*0061 % load 2-D and 3-D variable diagnostic output, annual mean data
0062 TRELAX = ncread('surfDiag.nc', 'TRELAX'); % (x, y, depth, time); depth=1 (surface-only)
0063 ETAN = ncread('surfDiag.nc', 'ETAN'); % (x, y, depth, time); depth=1 (surface-only)
0064 UVEL = ncread('dynDiag.nc', 'UVEL'); % (x, y, depth, time); x dim is Nx+1, includes eastern edge
0065 VVEL = ncread('dynDiag.nc', 'UVEL'); % (x, y, depth, time); y dim is Ny+1, includes northern edge
0066 THETA = ncread('dynDiag.nc', 'THETA'); % (x, y, depth, time)
0067
0068 % MITgcm time unit (dim='T') is seconds, so read new time series
0069 % for (dynStDiag) monthly time averages and convert into years.
0070 % Note that the MITgcm time array values correspond to time at the end
0071 % of the time-avg periods, i.e. subtract 1/24 to plot at mid-month.
0072 Tmon = ncread('dynStDiag.nc', 'T')/(86400*360) - 1/24;
0073 % and repeat to read time series for annual mean time output,
0074 % subtract 0.5 to plot at mid-year.
0075 Tann = ncread('surfDiag.nc', 'T')/(86400*360) - 0.5;
0076
0077
0078 %%%%%%%%%%% plot diagnostics %%%%%%%%%%%
0079 % No sizing of figures is specified below, will render at matlab default
0080 % to view properly, resize windows appropriately
94151a9b18 Jeff*0081
0082 % figure 4.6 - time series of global mean TRELAX and THETA by level
ce0d9af5ea Jeff*0083 %
0084 % plot THETA at MITgcm k levels 1,5,15 (SST, -305 m, -1705 m)
0085 klevs = [1, 5, 15];
0086
0087 figure
0088 subplot(221)
0089 plot(Tmon,squeeze(TRELAX_ave(1,1,:)),'b','LineWidth',[4])
0090 % To specify colors for specific lines, matlab 2019b adds a
0091 % new function colororder(); otherwise, plot the lines separately
0092 % (providing a unique color spec) or grab the plot handle
0093 % and change colors after plotting.
0094 grid on
0095 hold on
0096 title('Net Heat Flux into Ocean (TRELAX\_ave)','FontSize',18)
0097 set(gca,'Fontsize', 18)
0098 xlabel('Time (yrs)')
0099 ylabel('W/m^2')
0100 set(gca,'YLim', [-400 0])
0101 % Alternatively, a global mean area-weighted TRELAX (annual mean)
0102 % could be computed as follows, using HfacC(:,:,1), i.e. HfacC in
0103 % the surface layer, as a land-ocean mask.
0104 % First, compute total surface area of ocean points:
0105 total_ocn_area = sum(sum(rA .* HFacC(:,:,1)));
0106 % Broadcasting the arrays across the time axis makes this a bit awkward:
0107 TRELAX_ave_ann = squeeze(sum(sum(TRELAX(:,:,1,:) .* ...
0108 repmat(rA .* HFacC(:,:,1), [1 1 1 size(TRELAX,4)])))) / total_ocn_area;
0109 plot(Tann, TRELAX_ave_ann, 'm--', 'LineWidth', [2])
0110 subplot(223)
0111 plot(Tmon, squeeze(THETA_lv_ave(klevs,1,:)), 'LineWidth', [4])
0112 grid on
0113 title('Mean Potential Temp. by Level (THETA\_lv\_avg)', 'FontSize', 18)
0114 xlabel('Time (yrs)')
0115 ylabel('^oC')
0116 legendStrings = string(RC(klevs)) + ' m';
0117 legend(legendStrings)
0118 set(gca,'Fontsize', 16)
0119 subplot(224)
0120 plot(Tmon, squeeze(THETA_lv_std(klevs,:)), 'LineWidth', [4])
0121 grid on
0122 title('Std. Dev. Potential Temp. by Level (THETA\_lv\_std)', 'FontSize', 18)
0123 set(gca,'Fontsize', 16)
0124 xlabel('Time (yrs)')
0125 ylabel('^oC')
0126 legend(legendStrings)
0127
0128 % figure 4.7 - 2-D plot of TRELAX and contours of free surface height
0129 % (ETAN) at simulation end ( endTime = 3110400000. is t=100 yrs).
0130 %
0131 mask = HFacC;
0132 mask(mask==0) = NaN; % create a 3-D land mask, with land points set to NaN
0133 eta_masked = mask(:,:,1) .* ETAN(:,:,1,end);
0134 figure
0135 pcolor(Xp1(1:end-1), Yp1(1:end-1), TRELAX(:,:,1,end)')
0136 colorbar
0137 shading flat
0138 colormap(bluetored)
0139 set(gca,'CLim', [-250 250])
0140 hold on
0141 contour(X,Y,eta_masked',[-.6:.1:.6],'k')
0142 set(gca,'XLim', [0 60])
0143 set(gca,'YLim', [15 75])
0144 title({'Free surface height (contours, CI .1 m)';'and TRELAX (shading, W/m^2)'})
0145 xlabel('Longitude')
0146 ylabel('Latitude')
0147 % Note we have used routine pcolor here with Xp1, Yp1, which are the
0148 % locations of the lower left corners of grid cells
0149 % (here, length Nx+1,Ny+1 as they include the ending right and upper
0150 % locations of the grid, respectively). We don't pass pcolor the
0151 % +1 location, but matlab pcolor ignores the last row and column
0152 % when plotting, here conveniently land points.
0153 % Alternative one could plot shading using contourf with 'LineStyle'
0154 % set to 'none' with coordinates X and Y, the grid cell center
0155 % locations. Also note we mask the land values as NaN when contouring
0156 % the free surface height.
94151a9b18 Jeff*0157
ce0d9af5ea Jeff*0158 % figure 4.8 - barotropic streamfunction, plot at simulation end
0159 % (w/overlaid labeled contours)
0160 %
0161 ubt = squeeze(sum(UVEL .* permute(repmat(drF, [1 Nx+1 Ny size(UVEL,4)]), ...
0162 [2 3 1 4]), 3)); % depth-integrated u velocity
0163 psi = zeros(Nx+1, Ny+1, size(UVEL,4));
0164 psi(:,2:end,:) = cumsum(-ubt .* repmat(dyG, [1 1 size(UVEL,4)]), 2) ...
0165 / 1.e6; % compute streamfn in Sv (for each yr)
0166 % Note psi is computed and plotted at the grid cell corners and we
0167 % compute as dimensioned (Nx+1,Ny); as noted, UVEL contains an extra
0168 % data point in x, at the eastern edge. cumsum is done in y-direction.
0169 % We have a wall at southern boundary (i.e. no reentrant flow from
0170 % north), ergo psi(j=1) is accomplished by declaring psi
0171 % to be shape (Nx+1, Ny+1, size(UVEL,4)) where size(UVEL, 4) is the
0172 % size of the time dimension.
0173
0174 figure
0175 clabel(contourf(Xp1, Yp1, psi(:,:,end)', [-35:5:35], 'k'))
0176 colorbar
0177 colormap(blueyelred)
0178 set(gca, 'XLim', [0 60])
0179 set(gca, 'Ylim', [15 75])
0180 set(gca, 'Clim', [-35 35])
0181 title('Barotropic Streamfunction (Sv)')
0182 xlabel('Longitude')
0183 ylabel('Latitude')
94151a9b18 Jeff*0184
ce0d9af5ea Jeff*0185 % figure 4.9 - potential temp at depth and xz slice, at simulation end
0186 %
0187 % plot THETA at MITgcm k=4 (220 m depth) and j=15 (28.5N)
0188 klev = 4;
0189 jloc = 15;
0190
0191 theta_masked = mask(:,:,klev) .* THETA(:,:,klev,end);
0192 figure
0193 subplot(121)
0194 % again we use pcolor for the plan view and provide the
0195 % corner point locations XG,YG via Xp1,Yp1
0196 pcolor(Xp1(1:end-1), Yp1(1:end-1), THETA(:,:,klev,end)')
0197 shading flat
0198 hold on
0199 colorbar
0200 colormap(coolwarm)
0201 % but to overlay contours we provide the cell centers
0202 % and mask out boundary/land cells
0203 % also note we are passing contour 2-D coordinates XC,YC
0204 % so no transpose on THETA needed
0205 contour(XC, YC, theta_masked, [0:2:30], 'k')
0206 set(gca, 'XLim', [0 60])
0207 set(gca, 'Ylim', [15 75])
0208 set(gca, 'Clim', [0 30])
0209 title(['a) THETA ' num2str(RC(klev)) ' m Depth (^oC)'])
0210 xlabel('Longitude')
0211 ylabel('Latitude')
0212
0213 % For the xz slice, our limited vertical resolution makes for an ugly
0214 % pcolor plot, we'll shade using contour instead, providing the centers
0215 % of the vertical grid cells and cell centers in the x-dimension.
0216 % Also mask out land cells at the boundary, which results in slight
0217 % white space at domain edges.
0218 theta_masked = squeeze(mask(:,jloc,:).*THETA(:,jloc,:,end));
0219 subplot(122)
0220 contourf(X, RC, theta_masked', 0:2:30)
0221 colorbar
0222 title(['b) THETA at ' num2str(Y(jloc)) 'N (^oC)'])
0223 set(gca, 'CLim', [0 30])
0224 set(gca, 'XLim', [0 60])
0225 set(gca, 'YLim', [-1800 0])
0226 xlabel('Longitude')
0227 ylabel('Depth (m)')
0228 % One approach to avoiding the white space at boundaries/top/bottom
0229 % is to copy neighboring tracer values to land cells prior to contouring
0230 % (and don't mask), and augment a row of z=1 data at the ocean surface
0231 % and a row at the ocean bottom.
0232 % (see tutorial Southern Ocean Reentrant Channel for an example)
0233 % The white space occurs because contour uses data at cell centers, with
0234 % values masked/undefined beyond the cell centers toward boundaries.
0235
0236 % To instead plot using pcolor, pass location of vertical cell faces RF:
0237 %
0238 % RF=ncread('grid.nc', 'RF')
0239 % pcolor(Xp1(1:end-1), RF(1:end-1), squeeze(THETA(:,jloc,:,end))')
0240 % shading flat; colormap(coolwarm)
0241 %
0242 % but note that unlike python, matlab pcolor doesn't plot the bottom
0243 % layer of valid temperature data (the ending column)
