** Warning **
Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.
MITgcm/MITgcm/verification/tutorial_reentrant_channel/analysis/matlab_plots.m
Warning, /verification/tutorial_reentrant_channel/analysis/matlab_plots.m is written in an unsupported language. File is not indexed.
view on github raw file Latest commit 62748458 on 2021-08-11 00:15:22 UTC
62748458e3 Jeff* 0001 % matlab plotting code for Tutorial Southern Ocean Reentrant Channel
0002
0003 % read in additional colormaps for plots
0004 % assumes 'addpath ../../../utils/matlab/' typed from your run directory
0005 % which is also needed for other matlab functions such as 'rdmds'
0006
0007 bluered_colormaps
0008
0009 use_GM = 1; % set to zero if not using GM
7655a8a3af Jeff* 0010
0011 %%%%%%%%%%% load grid data
0012
62748458e3 Jeff* 0013 XC = rdmds('XC');
0014 YC = rdmds('YC');
0015 XG = rdmds('XG');
0016 YG = rdmds('YG');
0017 RC = squeeze(rdmds('RC')); % dimensioned (1,1,Nr) so squeeze to vector form
0018 RF = squeeze(rdmds('RF')); % dimensioned (1,1,Nr+1) so squeeze to vector form
0019 DRC = rdmds('DRC');
0020 DRF = rdmds('DRF');
0021 DXG = rdmds('DXG');
0022 DYG = rdmds('DYG');
0023 DXC = rdmds('DXC');
0024 DYC = rdmds('DYC');
0025 rAC = rdmds('rAC'); % gridcell area, cell centers, aka rA
0026 hFacC = rdmds('hFacC'); % partial cell hFac weightings: cell centers
0027 hFacW = rdmds('hFacW'); % partial cell hFac weightings: west face
0028 hFacS = rdmds('hFacS'); % partial cell hFac weightings: south face
7655a8a3af Jeff* 0029
62748458e3 Jeff* 0030 % number of gridcells in (x,y,z)
0031 Nx = size(XC,1);
0032 Ny = size(XC,2);
0033 Nr = length(RC);
7655a8a3af Jeff* 0034
62748458e3 Jeff* 0035 %create 2D XGp1, YGp1 to span full domain (0:1000 km, 0:2000 km)
0036 YGp1=[YG YG(:,end)+DYG(:,end)];
0037 YGp1=[YGp1; YGp1(end,:)];
0038 XGp1=[XG; XG(end,:)+DXG(end,:)];
0039 XGp1=[XGp1 XGp1(:,end)];
7655a8a3af Jeff* 0040
62748458e3 Jeff* 0041 % define layers as specified in data.layers
0042 Tlay = [-2.00, -1.75, -1.50, -1.25, -1.00, -0.75, -0.50, -0.25, 0.00, ...
0043 0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, ...
0044 2.50, 2.75, 3.00, 3.25, 3.50, 3.75, 4.00, 4.25, 4.50, ...
0045 5.00, 5.50, 6.00, 6.50, 7.00, 7.50, 8.00, 8.50, 9.00, ...
0046 9.50, 10.00];
7655a8a3af Jeff* 0047
0048
62748458e3 Jeff* 0049 %%%%%%%%%%% load diagnostics
0050 %
0051 % as specified in data.diagnostics:
0052 % state: THETA, UVEL, VVEL, WVEL, CONVADJ
0053 % 2D_diags: TRELAX, MXLDEPTH, ETAN
0054 % layDiag: LaVH1TH, LaHs1TH, LaVa1TH
0055 % GM_diags: GM_PsiX , GM_PsiY (if use_GM)
0056
0057 if Ny==40 % coarse resolution
0058 tt = 933120; % choose iter for time averaged yr 30 data
0059 state=rdmds('Diags/state', tt);
0060 diag2D=rdmds('Diags/2D_diags', tt);
0061 if use_GM
0062 gm_diags = rdmds('Diags/GM_diags', tt);
0063 end
0064 laydiag = rdmds('Diags/layDiag', tt);
7655a8a3af Jeff* 0065
62748458e3 Jeff* 0066 else %eddying run
0067
0068 % for eddy run, take a 5-year mean of yrs 26-30
0069 % instead of using single-yr mean
0070 ttrange = 3234816:124416:3732480;
0071
0072 state = zeros(Nx, Ny, Nr, 5);
0073 for tt = ttrange
0074 state=state+rdmds('Diags/state',tt);
0075 end
0076 state = state / length(ttrange);
0077
0078 diag2D = zeros(Nx,Ny,3);
0079 for tt = ttrange
0080 diag2D = diag2D + rdmds('Diags/2D_diags',tt);
0081 end
0082 diag2D = diag2D / length(ttrange);
0083
0084 laydiag = zeros(Nx,Ny,length(Tlay)-1,3);
0085 for tt = ttrange
0086 laydiag = laydiag + rdmds('Diags/layDiag',tt);
0087 end
0088 laydiag = laydiag / length(ttrange);
0089
0090 end
7655a8a3af Jeff* 0091
0092
0093 %%%%%%%%%%% plot diagnostics
0094
62748458e3 Jeff* 0095 % statistical diagnostics output
0096 %
0097 % first, parse and load data
0098 % assumes you have run extract_StD script, see tutorial
0099 [nIter, regList, time, stdiagout, listFlds, listK] = read_StD('STATDIAGS', 'dat', 'all_flds');
0100 % plot time series of TRELAX
0101 figure
0102 plot(time(:,2) / (86400*360), stdiagout(1,:,1,1,2), 'b', 'LineWidth', 4)
0103 grid on
0104 hold on
0105 plot(time(:,2) / (86400*360), stdiagout(1,:,1,2,2), 'm', 'LineWidth', 4)
0106 legend('Mean', 'Std Dev')
0107 title('Net Heat Flux into Ocean')
0108 xlabel('Time (yrs)')
0109 ylabel('W/m^2')
0110 set(gca,'FontSize', 14)
0111 % plot time series of THETA at surface, 270 m and 2580 m
0112 figure
0113 subplot(3,1,1) % MITgcm k=1 SST (region=1, field =1)
0114 plot(time(:,2) / (86400*360), stdiagout(2,:,1,1,1), 'c', 'LineWidth', 4)
0115 grid on
0116 title('Mean Temperature, Surface')
0117 xlabel('Time (yrs)')
0118 ylabel('^oC')
0119 set(gca,'FontSize', 14)
0120 subplot(3,1,2) % MITgcm k=17 270m depth (THETA)
0121 plot(time(:,2) / (86400*360), stdiagout(18,:,1,1,1), 'g', 'LineWidth', 4)
0122 grid on
0123 title('Mean Temperature, Depth 270 m')
0124 xlabel('Time (yrs)')
0125 ylabel('^oC')
0126 set(gca,'FontSize', 14)
0127 subplot(3,1,3) % MITgcm k=40 2577m depth (THETA)
0128 plot(time(:,2) / (86400*360), stdiagout(41,:,1,1,1), 'r', 'LineWidth', 4)
0129 grid on
0130 title('Mean Temperature, Depth 2580 m')
0131 xlabel('Time (yrs)')
0132 ylabel('^oC')
0133 set(gca,'FontSize', 14)
0134
7655a8a3af Jeff* 0135 % zonal mean temperture and mixed-layer depth
0136 %
62748458e3 Jeff* 0137 % weight temperature by partial cell factor hFacC to compute zonal mean
0138 zmT = squeeze(sum(state(:,:,:,1) .* hFacC .* rAC, 1)) ./ squeeze(sum(hFacC .* rAC, 1));
0139 % the next steps create a "padded" version of zmT which we will plot
0140 % to avoid undefined values outside cell centers at bottom and N,S edges
0141 % (which matlab will assume leave unshaded rather than extrapolate)
0142 % zmTpad is dimensioned w/extra value in y (north edge) and z (bottom)
0143 zmTpad = zeros(Ny+1, length(RF));
0144 zmTpad(1:Ny,1:end-1) = zmT;
0145 % for southern edge, simply copy the neighbor onto land point
0146 zmTpad(1,1:end-1) = zmT(2,:);
0147 % for northern edge, copy the temperature from cell immediately to south
0148 zmTpad(end,1:end-1) = zmT(end,:);
0149 % for bottom, copy the temperature from above onto new bottom point
0150 zmTpad(:,end) = zmTpad(:,end-1);
0151 figure
0152 contourf([YC(1,:) YC(1,end)+DYC(1,end)]/1000, [RC' RC(end)-DRC(end)], ...
0153 zmTpad', -2:.1:10, 'Linestyle', 'none')
0154 set(gca, 'CLim', [-2 10])
0155 colorbar
0156 hold on
0157 zmML = squeeze(mean(diag2D(:,:,2),1)); % take zonal mean of mixed layer depth
0158 % for a given (ocean) latitude band, all ML data are ocean w/same area
0159 % so no need for masking land or area weighting in computing mean
0160 if Ny==40 % coarse res
0161 plot(YC(1,2:end)/1000, -zmML(2:end), 'k', 'LineWidth', 3) % do not plot ML over land
0162 else % eddying resolution
0163 plot(YC(1,11:end)/1000, -zmML(11:end), 'k', 'LineWidth', 3)
0164 end
0165 % shade land as grey
0166 rectangle('position',[0 RF(end) 50 -RF(end)], 'Facecolor',[.7 .7 .7], 'LineStyle', 'none');
0167 set(gca, 'YLim', [-3982 0])
0168 set(gca, 'XLim', [0 2000])
0169 xlabel('y-coordinate (km)')
0170 ylabel('Depth (m)')
0171 set(gca, 'FontSize', 14)
0172 title({'Zonal Mean Temperature (^oC)', 'Mixed Layer Depth (m)'})
0173
0174 % convective adjustment index plot
0175 %
0176 % plot MITgcm k=4 (92 m depth)
0177 klev = 10;
0178 convadj = zeros(Nx+1, Ny+1);
0179 % note extra row and column of zeros for matlab pcolor (ignored in plot)
0180 convadj(1:end-1,1:end-1) = state(:,:,klev,5);
0181 figure
0182 pcolor(XGp1/1000, YGp1/1000, convadj)
0183 shading flat
0184 colormap(flipud(summer))
0185 colorbar
0186 set(gca, 'PlotBoxAspectRatio', [1 2 2])
0187 xlabel('x-coordinate (km)')
0188 ylabel('y-coordinate (km)')
0189 set(gca, 'FontSize', 14)
0190 rectangle('position', [0 0 1000 50], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none')
0191 title(['Convective Adj. Index (depth ' num2str(abs(round(RC(klev)))) ' m)'])
0192 % users manual figure used colormap cmocean_v2 'deep' instead of 'summer'
0193 % https://tos.org/oceanography/article/true-colors-of-oceanography-guidelines-for-effective-and-accurate-colormap
0194 % available for download at matlab file exchange
0195 % https://www.mathworks.com/matlabcentral/fileexchange/57773-cmocean-perceptually-uniform-colormaps?s_tid=srchtitle
7655a8a3af Jeff* 0196
0197 % barotropic streamfunction
0198 %
62748458e3 Jeff* 0199 % see tutorial Baroclinic Gyre analysis for additional details on computation
0200 % given our periodic EW configuration, repeating i=1 data at i=end+1
7655a8a3af Jeff* 0201 % first, compute depth-integrated u velocity weighted by partial cell factor hFacW
62748458e3 Jeff* 0202 ubt = squeeze(sum(state([1:end 1],:,:,3) .* hFacW([1:end 1],:,:) .* repmat(DRF, [Nx+1 Ny 1]), 3));
7655a8a3af Jeff* 0203 % next, need to include row of zeros at y=0, then cumsum -ubt.*DYG in y direction
62748458e3 Jeff* 0204 psi = [zeros(Nx+1,1) cumsum(-ubt .* DYG([1:end 1],:), 2)]/1.e6; % and convert to Sv
0205 figure
0206 contourf(XGp1'/1000, YGp1'/1000, psi', [-360:20:0])
0207 colorbar
0208 colormap('summer');
0209 % tutorial figure used colormap cmocean_v2 '-speed' instead of matlab 'summer'
0210 rectangle('position', [0 0 1000 50], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none');
0211 set(gca,'CLim',[-350 0])
0212 set(gca,'PlotBoxAspectRatio', [1 2 2])
0213 set(gca,'XLim', [0 1000])
0214 set(gca,'YLim', [0 2000])
0215 xlabel('x-coordinate (km)')
0216 ylabel('y-coordinate (km)')
0217 set(gca,'FontSize', 14)
0218 title('Barotropic Streamfunction (Sv)')
7655a8a3af Jeff* 0219
0220 % Eulerian MOC
0221 %
62748458e3 Jeff* 0222 % first, take zonal sum of v*dx*dz, adding a row of zeros at the bottom
0223 % (RF spans 0:-3982 m dimension 50)
0224 vzi = [ squeeze(sum(state(:,:,:,2) .* hFacS .* repmat(DXG, [1 1 length(RC)]) ...
0225 .* repmat(DRF, [Nx Ny]), 1)) zeros(Ny, 1) ];
0226 % next, best in general to do the cumsum bottom up
0227 % avoids issues if one uses a nonlinear free surface
0228 moc = -cumsum(vzi, 2, 'reverse')/1.e6; % and convert to Sv
0229 figure
0230 contourf(YG(1,:)/1000, RF, moc', [-5:.1:5], 'Linestyle', 'none')
0231 set(gca, 'Clim', [-3 3])
0232 hold on
0233 contour(YG(1,:)/1000, RF, moc', [-5:.5:5], 'k')
0234 set(gca, 'Clim', [-3 3])
0235 colorbar
0236 rectangle('position', [0 RF(end) 50 -RF(end)], 'Facecolor', [.7 .7 .7], 'LineStyle','none');
0237 xlabel('y-coordinate (km)')
0238 ylabel('Depth (m)')
0239 set(gca,'FontSize',14)
0240 colormap(bluetored)
0241 title('Eulerian MOC (Sv)');
7655a8a3af Jeff* 0242
0243 % Eulerian MOC plus Bolus velocity
0244 %
62748458e3 Jeff* 0245 % to add bolus velocity to Eulerian MOC (i.e., residual MOC), if using GM
0246 % compute bolus velocity from bolus streamfunctionn diagnostic
0247 % (using advective aka bolus form of GM)
0248 if use_GM
0249 psiy = zeros(Nx, Ny, Nr+1);
0250 % add a plane of zeros to psi at bottom of domain, we will take derivative in z
0251 psiy(:,:,1:end-1) = gm_diags(:,:,:,2);
0252 bolV = (psiy(:,:,2:end) - psiy(:,:,1:end-1)) ./ repmat(DRF, [Nx Ny]);
0253 % include bolV when computing vzi
7655a8a3af Jeff* 0254 % note bolV already has hFacS scaling factored into psiy (a transport diagnostic)
62748458e3 Jeff* 0255 vzi = [ squeeze(sum((state(:,:,:,2) .* hFacS + bolV) .* repmat(DXG, [1 1 Nr]) ...
0256 .* repmat(DRF, [Nx Ny]),1)) zeros(Ny, 1) ];
0257 moc = -cumsum(vzi, 2, 'reverse')/1.e6;
0258 figure
0259 contourf(YG(1,:)/1000, RF, moc', [-5:.1:5], 'Linestyle', 'none')
0260 set(gca, 'Clim', [-3 3])
0261 hold on;
0262 contour(YG(1,:)/1000, RF, moc', [-5:.5:5], 'k')
0263 set(gca, 'Clim', [-3 3])
0264 colorbar
0265 rectangle('position', [0 RF(end) 50 -RF(end)], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none')
0266 xlabel('y-coordinate (km)')
0267 ylabel('Depth (m)')
0268 set(gca,'FontSize',14)
0269 colormap(bluetored)
0270 title('Eulerian plus Bolus MOC (Sv)')
0271 end
7655a8a3af Jeff* 0272
0273 % residual MOC from layers, in density space
0274 %
0275 % first, using diagnostic LaVH1TH -- the meridional mass transport in layers
0276 % do a cumsum of the zonally integrated transport and multiply by DXG
0277 % note cumsum done forward order in matlab (1:37) but this is bottom-up in the ocean
0278 % e.g layer 1 is between -2 deg C and -1.75 deg C
62748458e3 Jeff* 0279 vti = [ zeros(Ny, 1)'; -repmat(DXG(1,:), [length(Tlay)-1 1]) ...
0280 .* cumsum(squeeze(sum(laydiag(:,:,:,1), 1)), 2)'/1.e6 ];
0281 % due to error in calculation, at the ocean top the value
0282 % is not precisely zero; force small residual -> 0
0283 vti(abs(vti)<.005) = 0;
7655a8a3af Jeff* 0284 % next, plot in density (temperature) space
62748458e3 Jeff* 0285 figure
0286 contourf(YG(1,:)/1000, Tlay, vti, [-5:.01:5], 'Linestyle', 'none')
0287 set(gca,'Clim', [-3 3])
0288 hold on
0289 contour(YG(1,:)/1000, Tlay, vti, [-3:.5:3], 'k')
0290 colorbar
7655a8a3af Jeff* 0291 % finally, plot bounds of max and min SST
62748458e3 Jeff* 0292 sstmax = max(state(:,:,1,1));
0293 sstmin = min(state(:,:,1,1)); % max, min in each latitude band
0294 if Ny==40
0295 plot(YC(1,2:end)/1000, sstmin(2:end), 'g--', 'LineWidth', 2)
0296 plot(YC(1,2:end)/1000, sstmax(2:end), 'g--', 'LineWidth',2)
0297 rectangle('position', [0 -2 50 sstmax(2)+2], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none')
0298 end
0299 % for eddy run, need to dump data more frequently to obtain reasonable
0300 % estimate for SST max, min (e.g. daily)
0301 xlabel('y-coordinate (km)')
0302 ylabel('Temperature (^oC)')
0303 set(gca,'FontSize',[14])
0304 colormap(bluetored)
0305 title('Residual Overturning Circulation (Sv)')
7655a8a3af Jeff* 0306
0307 % residual MOC from layers, converted back into depth coordinates
0308 %
0309 % to compute the z-locations of the layer interfaces, use diagnostic LaHs1TH (layer depths)
62748458e3 Jeff* 0310 % do a cumsum from the ocean surface downward
0311 % (thus we flip the laydiag data in z-axis so that k=1 is surface)
7655a8a3af Jeff* 0312 % and then take a zonal mean of these z-locations (obviously, need to be
0313 % careful about this if the layer depths vary radically in a zonal band!)
62748458e3 Jeff* 0314 csum_th = -[ zeros(Ny, 1) squeeze(mean(cumsum(laydiag(:,:,end:-1:1,2), 3), 1)) ];
0315 % then cumsum the layer transports from ocean surface downward,
0316 % multiply by DXG, and zonally integrate
0317 csum_tr = [ zeros(Ny, 1) squeeze(sum(cumsum(laydiag(:,:,end:-1:1,1), 3) ...
0318 .* repmat(DXG, [1 1 length(Tlay)-1]), 1)) ]/1.e6;
0319 % then this can be plotted specifying 2D arrays for x-axis (Y-coordinate)
0320 % and y-axis (layer z-locations)
0321 Y = repmat(YG(1,:), [length(Tlay) 1])';
0322 figure
0323 rectangle('position', [0 -4000 2000 520.5], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none')
0324 hold on
0325 contourf(Y/1000, csum_th, csum_tr, [-5:.01:5], 'Linestyle', 'none')
0326 colorbar
0327 set(gca, 'CLim', [-3 3])
0328 % plot zonal mean temp contours, zmT computed above
0329 contour(YC(1,:)/1000, RC, zmT', -2:1:10, 'k');
0330 xlabel('y-coordinate (km)')
0331 ylabel('Depth (m)')
0332 set(gca,'FontSize',14)
0333 colormap(bluetored);
0334 rectangle('position', [0 -4000 50 4000], 'Facecolor', [.7 .7 .7], 'LineStyle', 'none');
0335 title('Residual MOC converted to depth space (Sv)')
0336 set(gca,'YLim', [-4000 0])
