do_volume_change2.m 4.84 KB
Newer Older
 Pierre Cazenave committed Jun 20, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 function [summedVolume]=do_volume_change2(plotOPTS,FVCOM,startIdx) % Calculate the volume exposed to a given change in pH. m_mappath; [nx,nz,ttot]=size(FVCOM.(plotOPTS.var_plot)); % Check we have some of the required fields. if ~isfield(FVCOM,plotOPTS.var_plot) error('Need %s input to calculate change in %s.',plotOPTS.var_plot,plotOPTS.var_plot) end % Get the time step (in seconds) dt=round((plotOPTS.Time_record(2)-plotOPTS.Time_record(1))*24*60*60); % Sigma layer fraction (nz*nz) sigThickness=roundn(abs(diff(FVCOM.siglev,1,2)),-5); % roundn to even values out. % Total depth (nz*ttot) totalDepth=repmat(FVCOM.h,1,size(FVCOM.zeta,2))+FVCOM.zeta; % Volume for every element for each time step cellVolume=nan(nx,nz,ttot); for ii=1:nz for jj=1:ttot % I'm pretty certain this needs to integrate for the time between % time steps, otherwise halving the output time step would halve % the exposed volume, which can't be right... cellVolume(:,ii,jj)=totalDepth(:,jj).*sigThickness(:,ii).*FVCOM.art1*dt; end end dataToUse=single(FVCOM.(plotOPTS.var_plot)); % Old way. % % Calculate the difference between the current step and the previous step. % phChange=diff(squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,:)),[],2); % % Nope, get the background condition and compare against that. Depth % average if necessary. % if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average % bgPH=squeeze(mean(FVCOM.(plotOPTS.var_plot)(:,:,startIdx),2)); % else % bgPH=squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,startIdx)); % end % % Are we depth averaging? % if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average % % Now, for each successive time step, calculate the difference between the % % current time step and the background level. % phMean=squeeze(mean(FVCOM.(plotOPTS.var_plot),2)); % phDiff=phMean-repmat(bgPH,1,size(phMean,2)); % else % phDiff=squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,:))-repmat(bgPH,1,size(FVCOM.(plotOPTS.var_plot),3)); % end % New, simple way. Just compare against the background. % Get the background condition. Depth average if necessary. if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average bgPH=squeeze(mean(dataToUse(:,:,startIdx),2)); else bgPH=squeeze(dataToUse(:,plotOPTS.nz_plot,startIdx)); end % Are we depth averaging? if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average phData=squeeze(mean(dataToUse,2)); else phData=squeeze(dataToUse(:,plotOPTS.nz_plot,:)); end % Get the differences phDiff=phData-repmat(bgPH,1,size(phData,2)); % For positions if we want to plot them [X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on'); foundXY=cell(1,length(plotOPTS.nz_plot)); countSites=zeros(nx,1); % Find all locations which experience a change in excess of the threshold. for tt=1:ttot % Find drops greater than the threshold idx=find(phDiff(:,tt)