** 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.

Last-Modified: Fri, 31 Oct 2024 05:11:45 GMT Content-Type: text/html; charset=utf-8 MITgcm/MITgcm/verification/tutorial_reentrant_channel/analysis/matlab_plots.m
Back to home page

MITgcm

 
 

    


Warning, /verification/tutorial_reentrant_channel/analysis/matlab_plots.m is written in an unsupported language. File is not indexed.

view on githubraw 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])