Update with the CO2 modelling scripts and functions

fvcom_postproc/plots/ph/fixmatlabfigures.sh

+#!/bin/bash
+
+# Use ImageMagick to convert ugly MATLAB pdfs to beautiful pngs.
+
+inFiles=("$@")
+
+for ((i=0; i<${#inFiles[@]}; i++)); do
+    echo -n "Fixing ${inFiles[i]}... "
+    convert \
+        -trim \
+        -density 600x600 \
+        "${inFiles[i]}" \
+        "${inFiles[i]%.*}".png
+    echo "done."
+done

fvcom_postproc/plots/volumes.csv

+filename,"volume (m3) -0.2","volume (m3) -0.5","volume (m3) -1.0"
+co2_S5_high_run_fvcom_inputV5_low_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_low_run_fvcom_inputV5_low_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_pipe_run_fvcom_inputV5_low_flow_0001.nc,37049053187475.59,9401069500260.20,4038635613015.17,
+co2_S5_high_run_fvcom_inputV5_high_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_low_run_fvcom_inputV5_high_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_pipe_run_fvcom_inputV5_high_flow_0001.nc,55180178565003.48,16709666695708.14,8426030926795.49,
+co2_S5_high_run_fvcom_inputV5_mid_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_low_run_fvcom_inputV5_mid_flow_0001.nc,0.00,0.00,0.00,
+co2_S5_pipe_run_fvcom_inputV5_mid_flow_0001.nc,15221674579876.93,6047205278629.37,1666350432064.81,
+co2_S7_low_run_fvcom_inputV7_low_flow_0001.nc,0.00,0.00,0.00,
+co2_S7_high_run_fvcom_inputV7_high_flow_0001.nc,0.00,0.00,0.00,
+co2_S7_high_run_fvcom_inputV7_mid_flow_0001.nc,0.00,0.00,0.00,
+co2_S7_low_run_fvcom_inputV7_mid_flow_0001.nc,0.00,0.00,0.00,

fvcom_postproc/plots/volumes_temp.csv

+filename,"volume (m3) -0.2","volume (m3) -0.5","volume (m3) -1.0"

fvcom_postproc/surface_plots/get_velocity.m

+function [velMin,velMax,velMed,velMean,velStd]=get_velocity(u,v)
+% Get velocity stats from u and v components.
+% 
+% 
+
+vel = sqrt(u(:).^2+v(:).^2);
+velMin = min(vel);
+velMax = max(vel);
+velMed = median(vel);
+velMean = mean(vel);
+velStd = std(vel);

fvcom_postproc/surface_plots/scratch_figures.m

+close all
+
+plot(allIdx,testChange)
+hold on
+
+plot(allIdx(negIdx),negIdx(negIdx==1)*0,'.')
+plot(allIdx(negIdx),negIdx(negIdx==1)*0,'r.')
+
+plot(allIdx(idxStart_dn),negIdx(idxStart_dn)*-0.001,'g.')
+plot(allIdx(idxEnd_dn),negIdx(idxEnd_dn)*0.001,'k.')
+
+% plot(allIdx(idxStart_tc),negIdx(idxStart_tc)*-0.002,'gx')
+% plot(allIdx(idxEnd_tc),negIdx(idxEnd_tc)*0.002,'kx')
+
+% ylim([-5e-3 5e-3])
+axis([150 200 -5e-3 5e-3])
+
+%%
+
+for i=1:length(allInc)
+    figure(1)
+    clf
+    plot(testTime,squeeze(phChange(allInc(i),1,1:end-1)))
+    hold on
+    for j=1:size(longOnesIdx{allInc(i)},1)
+        plot(testTime(longOnesIdx{allInc(i)}(j,:)),zeros(2,1),'g-x')
+    end
+    text(min(testTime),0,num2str(allInc(i)))
+    ylim([-0.0003 0.0003])
+    pause(0.1)
+end
+
+%%
+
+for i=1:length(allInc)
+    figure(1)
+    clf
+    plot(testTime,squeeze(phChange(allInc(1),1,1:end-1)))
+    hold on
+    for j=1:size(longOnesIdx{1},1)
+        plot(testTime(longOnesIdx{1}(j,:)),zeros(2,1),'g-x')
+    end
+    text(min(testTime),0,num2str(allInc(1)))
+    ylim([-0.0003 0.0003])
+    pause(0.1)
+end
+
+%%
+
+figure(1)
+clf
+% plot(testTime,testData(testIdx,2:end))
+plot(testTime,testData(2:end))
+hold on
+for j=1:size(longOnesIdx,1)
+    plot(testTime(longOnesIdx(j,:)),zeros(2,1),'g-x')
+end

fvcom_postproc/surface_plots/total_CO2.m

+% Get the total CO2 for a day of model.
+close all
+
+leakidx=1316;
+dt=3600;
+
+[xt,zt,ttot]=size(FVCOM.DYE);
+
+nt=24; % Let's do a day.
+nBuff=12; % half day buffer
+dtStartIdx=(3600*nBuff)/dt; % in case we don't have hourly sampled output
+dtEndIdx=((3600*nt)/dt)+dtStartIdx;
+
+% 3D array of the layer thicknesses
+thickFactor=repmat(abs(diff(FVCOM.siglev,1,2)),[1,1,ttot]);
+% 3D array of the total depths at each element (tide + still water depth)
+totalDepth=permute(repmat(FVCOM.zeta+repmat(FVCOM.h,1,ttot),[1,1,zt]),[1,3,2]);
+% Now we have two 3D arrays, we just need to multiply totalDepth by
+% thickFactor to get layer thickness in metres
+layerThickness=thickFactor.*totalDepth;
+
+% Now, integrate the CO2 by depth across a day's worth of time steps.

fvcom_postproc/surface_plots/total_volume.m

+% Daily CO2 input for bottom cell
+
+close all
+clc
+
+leakidx=1316;
+baseidx=1;
+dt=3600;
+
+% Get bottom element input for 24 hours (skipping the buffer at the
+% beginning of half a day or so)
+nt=24; % Let's do a day.
+nBuff=12; % half day buffer
+dtStartIdx=(3600*nBuff)/dt; % in case we don't have hourly sampled output
+dtEndIdx=((3600*nt)/dt)+dtStartIdx;
+
+% Should just be the input amount at the bottom cell times the time step
+% times the number of time steps in a day.
+fprintf('Daily input (pica): %g\n',sum(squeeze(FVCOM.DYE(leakidx,baseidx,dtStartIdx:dtEndIdx))*dt))
+
+
+% My version: ((c*dt*n)/v)*86400
+fprintf('Daily input (pica): %g\n',(sum(squeeze(FVCOM.DYE(leakidx,baseidx,dtStartIdx:dtEndIdx)))*dt/(FVCOM.art1(leakidx)*abs(diff(FVCOM.siglev(leakidx,baseidx:baseidx+1)))))*86400)
+% Riqui's version: c*dt*n*v
+fprintf('Daily input (rito): %g\n',sum(squeeze(FVCOM.DYE(leakidx,baseidx,dtStartIdx:dtEndIdx)))*dt*FVCOM.art1(leakidx)*abs(diff(FVCOM.siglev(leakidx,baseidx:baseidx+1))))
+
+%% A new dawn...
+
+home
+
+% Modelled
+modelledInput=zeros(dtEndIdx-dtStartIdx+1,1);
+dtIter=dtStartIdx:dtEndIdx;
+for i=1:length(dtIter)
+    cellVolume=FVCOM.art1(leakidx)*((FVCOM.h(leakidx)+FVCOM.zeta(leakidx,dtIter(i)))*abs(diff(FVCOM.siglev(leakidx,baseidx:baseidx+1))));
+    if i>1
+        modelledInput(i)=modelledInput(i-1)+(((FVCOM.DYE(leakidx,baseidx,dtIter(i))*dt)/cellVolume));
+    end
+end
+TCO2=cumsum((squeeze(FVCOM.DYE(leakidx,baseidx,dtStartIdx:dtEndIdx))*dt)/cellVolume);
+
+% Theoretical
+theoreticalInput=zeros(dtEndIdx-dtStartIdx+1,1);
+for i=1:length(dtIter)
+    % Use the real cell volumes
+    cellVolume=FVCOM.art1(leakidx)*((FVCOM.h(leakidx)+FVCOM.zeta(leakidx,dtIter(i)))*abs(diff(FVCOM.siglev(leakidx,baseidx:baseidx+1))));
+    if i>1
+        theoreticalInput(i)=theoreticalInput(i-1)+(((500*dt)/cellVolume));
+    end
+end
+
+fprintf('Daily input (pica2): %g\n',modelledInput(end))
+
+fprintf('Daily input (theory): %g\n',theoreticalInput(end))
+fprintf('Daily input (theory2): %g\n',TCO2(end))
+
+
+%%
+
+% Compare the theoretical input (i.e. the input rate by the number of time
+% steps) against the modelled input at the leak point.
+close all
+
+clear modelledTCO2 theoreticalTCO2
+
+leakidx=1316;
+baseidx=1;
+dt=3600;
+inputCO2=500;
+
+colours=hsv(length(baseidx));
+
+nt=24; % Let's do a day.
+nBuff=12; % half day buffer
+dtStartIdx=(3600*nBuff)/dt; % in case we don't have hourly sampled output
+dtEndIdx=((3600*nt)/dt)+dtStartIdx;
+
+% Get the element area
+elemArea=FVCOM.art1(leakidx);
+
+% Get the modelled values
+for i=1:length(baseidx);
+    % Get the height of the base element
+    elemHeight=(FVCOM.zeta(leakidx,dtStartIdx:dtEndIdx)+FVCOM.h(leakidx)).*abs(diff(FVCOM.siglev(leakidx,baseidx(i):baseidx(i)+1)));
+    % Get the modelled CO2, and adjust for volume of base element
+    modelledTCO2(:,i)=(cumsum(squeeze(FVCOM.DYE(leakidx,baseidx(i),dtStartIdx:dtEndIdx)))./(elemArea.*elemHeight'))*dt;
+end
+
+% And the theoretical values (with the modelled tide effect i.e. change in
+% volume).
+theoreticalTCO2=((1:dtEndIdx-dtStartIdx+1)*500./(elemArea.*elemHeight))*dt;
+
+figure(1)
+hold on
+for i=1:length(baseidx);
+    plot(Time_record(dtStartIdx:dtEndIdx),modelledTCO2,'-x','LineWidth',3)
+end
+plot(Time_record(dtStartIdx:dtEndIdx),theoreticalTCO2,'r--x','LineWidth',2)
+axis('tight')
+ylabel('DYE')
+xlabel('Time (days)')
+legend('Modelled','Theoretical','Location','NorthWest')
+legend('BoxOff')
\ No newline at end of file

fvcom_postproc/surface_plots/volume_calculation.m

+% Time to sort out this volume calculation. 
+
+% For a given concentration going in (say 500 mmol/s), then the total
+% amount enetered into the system will be less because it'll be spread over
+% the volume of the TCE. So, 
\ No newline at end of file

fvcom_postproc/transect_plots/example_transect_plot.m

+% Sample script to extract and generate transect plots of tracer variables 
+% Modify to suit your requirements
+%
+% DESCRIPTION:
+%    Extracts data from FVCOM and generates vertical sliced contour plots 
+%
+% INPUT: 
+%   
+%   
+%           
+% OUTPUT:
+%   
+% EXAMPLE USAGE
+%    See scripts and modify at will 
+%
+% Author(s):  
+%    Ricardo Torres Plymouth Marine Laboratory
+%
+% Revision history
+%   
+%==============================================================================
+
+clear
+close all
+clc
+
+addpath ../../utilities
+%%-----------------------------------------------------------
+% set directories default values
+%
+FVCOM_root = '';
+FVCOM_data_dir='/home_nfs/rito/models/FVCOM/runCO2_leak/output/';
+FVCOM_mat_dir='../mat/';
+FVCOM_plot_dir='../plots/';
+FVCOM_ ='./';
+%
+time_offset = 678942; % from FVCOM time to matlab time
+%%
+%------------------------------------------------------------------------------
+% set casename specifics here
+%------------------------------------------------------------------------------
+%
+casename='co2_S5';
+experiment_name='slowleak';
+base_year = datenum(2006,1,1,0,0,0);
+files_FVCOM ='co2_S5.1.2.1_0002.nc';%
+date_range={'30/01/06 00:00:00','01/02/06 23:00:00'}; %'19/03/06 23:00:00' -1 if all available data wanted 
+plotOPTS.fig_name = [casename '_' experiment_name];
+var_2_xtractFVCOM = {'Itime','Itime2','xc','yc','h','siglay','siglev','zeta','TCO2','ph','u','v','temp','density'};
+dt=1/24; % time step of output in days.
+% Time record to extract
+%
+STnum = datenum(date_range{1},'dd/mm/yy HH:MM:SS');
+ENDnum = datenum(date_range{2},'dd/mm/yy HH:MM:SS');
+Time_record=STnum:dt:ENDnum;
+CD=pwd;
+%%
+%------------------------------------------------------------------------------
+% Specify What indices to read from FVCOM file
+%------------------------------------------------------------------------------
+siglev_idx=-1; % to extract all water levels from netcdf file
+siglay_idx=-1; % to extract all water levels from netcdf file
+
+%------------------------------------------------------------------------------
+load(fullfile(FVCOM_mat_dir, [casename, 'mesh']));
+plotOPTS.mesh=mesh;
+%------------------------------------------------------------------------------
+% select transect nodes with GUI
+%------------------------------------------------------------------------------
+trn_nodes= transect_nodes_screen(mesh,'nodes')
+node_idx=trn_nodes.idx;
+nele_idx=-1;
+%------------------------------------------------------------------------------
+% calculate distance along transect
+% ------------------------------------------------------------------------
+trn_dist=[trn_nodes.x(1)-trn_nodes.x(:),...
+    trn_nodes.y(1)-trn_nodes.y(:)];
+trn_dis=cumsum(abs(complex(trn_dist(:,1),trn_dist(:,2))));
+%
+%%
+%------------------------------------------------------------------------------
+% Specify variable and conditions of plot
+%------------------------------------------------------------------------------
+plotOPTS.var_plot = 'ph';
+plotOPTS.clims=[0.026418 0.026419];
+plotOPTS.do_mesh=0; % don't display mesh
+plotOPTS.trn_dis=trn_dis;
+plotOPTS.trn_nodes=trn_nodes;
+[plotOPTS.range_lat ,plotOPTS.range_lon] = deal([50.1 50.4],[ -4.5 -3.85]);
+plotOPTS.coastline_file=[FVCOM_mat_dir 'tamar3_0coast.mat'];
+plotOPTS.zone=30;
+plotOPTS.ell='grs80';
+plotOPTS.do_mesh=0; % don't display mesh
+
+%% ------------------------------------------------------------------------
+% 
+%------------------------------------------------------------------------------
+% read FVCOM data
+%------------------------------------------------------------------------------
+% 
+FVCOM_data=read_netCDF_FVCOM('time',date_range,'data_dir',FVCOM_data_dir ,...
+    'file_netcdf',files_FVCOM,'node_idx',node_idx,'nele_idx',nele_idx,'siglev_idx',siglev_idx,...
+    'siglay_idx',siglay_idx,'varnames',var_2_xtractFVCOM)
+% put variables into a strcuture variable
+for vv=1:length(var_2_xtractFVCOM)
+    FVCOM.(var_2_xtractFVCOM{vv}) = FVCOM_data{vv};
+end
+
+cd(CD)
+plotOPTS.figure=1;
+[Plots]=do_transect_plot(plotOPTS,FVCOM)
+
+
+

fvcom_postproc/transect_plots/get_transects.m

+% Get transects for the coarse and fine grids for the pH profiles
+
+% Load the results and the mesh
+
+clear
+close all
+clc
+
+% Redo the transect or load a saved one?
+reselecttransect='no';
+
+addpath ../../utilities
+addpath /users/modellers/rito/Models/fvcom-toolbox/utilities/
+%-----------------------------------------------------------
+% set directories default values
+%
+FVCOM_root = '';
+% FVCOM_data_dir='/data/medusa/rito/models/FVCOM/runCO2_leak/output/';
+FVCOM_data_dir='/data/medusa/pica/models/FVCOM/runCO2_leak/output/scenarios/';
+% FVCOM_data_dir='/tmp/data/20days/';
+% FVCOM_data_dir='/data/medusa/pica/models/FVCOM/runCO2_leak/output/rate_ranges/';
+% FVCOM_data_dir='/data/medusa/pica/models/FVCOM/runCO2_leak/output/sponge_tests/';
+FVCOM_mat_dir='/tmp/fvcom_co2/mat/';
+plotOPTS.FVCOM_plot_dir='/tmp/fvcom_co2/plots/';
+FVCOM_ ='./';
+%
+time_offset = 678942; % from FVCOM time to matlab time
+%
+%------------------------------------------------------------------------------
+% set casename specifics here
+%------------------------------------------------------------------------------
+[plotOPTS.range_lat ,plotOPTS.range_lon] = deal([50.2 50.4],[ -4.34 -3.91]);
+base_year = datenum(2006,1,1,0,0,0);
+
+% Use two results only -- doesn't matter what co2 scenario, we're only
+% interested in nodes and z
+allFiles={
+    'co2_S5_high_run_fvcom_inputV5_low_flow_0001.nc',...   % Coarse, fast leak
+    'co2_S7_low_run_fvcom_inputV7_low_flow_0001.nc'...    % Fine, slow leak
+};
+
+for fileIdx=1%:size(allFiles,2)
+    
+    files_FVCOM = allFiles{fileIdx};
+
+    if ~isempty(strfind(files_FVCOM,'S5')) || ~isempty(strfind(files_FVCOM,'V5'))
+        saveFile=fullfile(FVCOM_mat_dir,'S5_east-west_transect.mat');
+    elseif ~isempty(strfind(files_FVCOM,'S7')) || ~isempty(strfind(files_FVCOM,'V7'))
+        saveFile=fullfile(FVCOM_mat_dir,'S7_east-west_transect.mat');
+    else
+        error('Unknown grid resolution type. Can''t set output file name.')
+    end
+
+    
+    if exist(fullfile(FVCOM_data_dir,files_FVCOM), 'file')~=2
+        error('File not found %s',fullfile(FVCOM_data_dir,files_FVCOM))
+    end
+    
+    % Scenarios excluding warm up
+    date_range={'05/01/06 22:00:00','19/01/06 00:00:00'}; %'19/03/06 23:00:00' -1 if all available data wanted 
+
+    casename=files_FVCOM(1:6); % 'co2_S?'
+
+    % Output figure base name
+    % plotOPTS.fig_name = [casename '_' experiment_name];
+    plotOPTS.fig_name = files_FVCOM(1:end-3);
+
+%     var_2_xtractFVCOM = {'Itime','Itime2','xc','yc','h','siglay','siglev','zeta','TCO2','ph','art1','DYE'};
+    var_2_xtractFVCOM = {'Itime','Itime2','xc','yc','h','siglay','siglev','zeta','TCO2','ph','art1','DYE','u','v'};
+    dt=1/24; % time step of output in days.
+
+    % Time record to extract
+    %
+    STnum = datenum(date_range{1},'dd/mm/yy HH:MM:SS');
+    ENDnum = datenum(date_range{2},'dd/mm/yy HH:MM:SS');
+    Time_record=STnum:dt:ENDnum;
+    CD=pwd;
+    %
+    %------------------------------------------------------------------------------
+    % Specify what indices to read from FVCOM file
+    %------------------------------------------------------------------------------
+    siglev_idx=-1; % to extract all water levels from netcdf file
+    siglay_idx=-1; % to extract all water levels from netcdf file
+
+    %
+    %------------------------------------------------------------------------------
+    % load mesh information for current casename
+    %------------------------------------------------------------------------------
+    load(fullfile(FVCOM_mat_dir, [casename, 'mesh']));
+    plotOPTS.mesh=mesh;
+
+    %% Select the transect
+    if strcmpi(reselecttransect,'yes')
+        trn_nodes=transect_nodes_screen(mesh,'nodes')
+    else
+        trn_nodes=load(saveFile);
+        trn_nodes=trn_nodes.trn_nodes;
+    end
+    node_idx=trn_nodes.idx;
+    nele_idx=-1;
+    
+    %
+    %------------------------------------------------------------------------------
+    % Specify variable and conditions of plot
+    %------------------------------------------------------------------------------
+    plotOPTS.coastline_file=[FVCOM_mat_dir 'tamar3_0coast.mat'];
+    plotOPTS.zone=30;
+    plotOPTS.Time_record=Time_record;
+    plotOPTS.ell='grs80';
+    plotOPTS.do_mesh=0; % don't display mesh
+    plotOPTS.var_plot = 'ph';
+    plotOPTS.trn_nodes=trn_nodes;
+
+    % Get distance along line
+    xNorm=trn_nodes.x-min(trn_nodes.x);
+    yNorm=trn_nodes.y-min(trn_nodes.y);
+    xDist=sqrt(xNorm.^2+yNorm.^2);
+    plotOPTS.trn_dis=xDist;    
+    
+    %------------------------------------------------------------------------------
+    % 
+    %------------------------------------------------------------------------------
+    % read FVCOM data
+    %------------------------------------------------------------------------------
+    % 
+    FVCOM_data=read_netCDF_FVCOM('time',date_range,'data_dir',FVCOM_data_dir ,...
+        'file_netcdf',files_FVCOM,'node_idx',node_idx,'nele_idx',nele_idx,'siglev_idx',siglev_idx,...
+        'siglay_idx',siglay_idx,'varnames',var_2_xtractFVCOM);
+    % put variables into a strcuture variable
+    for vv=1:length(var_2_xtractFVCOM)
+        FVCOM.(var_2_xtractFVCOM{vv}) = FVCOM_data{vv};
+    end
+
+
+    
+    % Check the position of the transect
+%     cd(CD)
+%     plotOPTS.figure=1;
+%     plotOPTS.trn_nodes=trn_nodes;
+%     [Plots]=do_transect_plot(plotOPTS,FVCOM);
+
+    % Check the depths we've extracted
+    figure
+    plot(xDist,-FVCOM.h)
+    
+    % Save the transect information to the save file
+    if strcmpi(reselecttransect,'yes')
+        save(saveFile,'trn_nodes')
+    end
+    
+    %% Save a figure of the transect location
+    
+    % Re-read the results
+    %------------------------------------------------------------------------------
+    % 
+    %------------------------------------------------------------------------------
+    % read FVCOM data
+    %------------------------------------------------------------------------------
+    % 
+    node_idx_all=-1; % to extract all nodes from netcdf file
+    nele_idx_all=-1; % to extract all elements from netcdf file
+
+    FVCOM_data=read_netCDF_FVCOM('time',date_range,'data_dir',FVCOM_data_dir ,...
+        'file_netcdf',files_FVCOM,'node_idx',node_idx_all,'nele_idx',nele_idx_all,'siglev_idx',siglev_idx,...
+        'siglay_idx',siglay_idx,'varnames',var_2_xtractFVCOM);
+    % put variables into a strcuture variable
+    for vv=1:length(var_2_xtractFVCOM)
+        FVCOM.(var_2_xtractFVCOM{vv}) = FVCOM_data{vv};
+    end
+
+    plotOPTS.clims=[-30 0];
+    m_mappath;
+    figure(100); clf
+    m_proj('UTM','lon',[plotOPTS.range_lon],'lat',[plotOPTS.range_lat],'zon',plotOPTS.zone,'ell','grs80')
+    m_grid('box','fancy')
+    m_usercoast(plotOPTS.coastline_file,'Color','k','LineWidth',3);
+    [X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
+    hold on
+    patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
+        'Cdata',-FVCOM.h,'edgecolor','interp','facecolor','interp');
+    caxis(plotOPTS.clims)
+    ch=colorbar;
+    colormap(jet)
+    ylabel(ch,'Depth (m)')
+    plot(X(node_idx), Y(node_idx),'-w.','LineWidth',2,'MarkerSize',5)
+%     plot(sort(X(node_idx)), sort(Y(node_idx),1,'descend'),'-w.','LineWidth',2,'MarkerSize',5)
+    
+    fprintf('Saving transect location figure... ')
+    set(findobj(gcf,'Type','text'),'FontSize',10)
+    %set(gcf,'PaperSize',fliplr(get(gcf,'PaperSize'))) 
+    set(gcf,'PaperPositionMode','auto');
+    set(gcf,'renderer','painters'); % for vector output in pdfs
+    print(gcf,'-dpdf','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/pdf/',casename,'_transect_location.pdf']); % pdf
+    %print(gcf,'-dpng','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/png/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_change_',num2str(plotOPTS.save_intervals(aa)),'.png']); % png
+    fprintf('done.\n')
+
+    
+    
+end
+

utilities/do_ph_change_plot.m

+function [Plots]=do_ph_change_plot(plotOPTS,FVCOM,startIdx)
+% Calculate the change in the given parameter (plotOPTS.var_plot) and plot
+% accordingly.
+m_mappath;
+
+% For testing
+% startIdx=48; % Leak start
+% layerIdx=10; % Bottom layer
+% For testing
+
+% 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 relevant time intervals
+if isempty(plotOPTS.save_intervals)
+    plotOPTS.save_intervals=1:length(plotOPTS.Time_record);
+    dontSave=1;
+else
+    dontSave=0;
+end
+
+% Build a colour palette which matches Jerry's ranges.
+% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
+nColours=200;
+nColourIn=[1,nColours*0.15,nColours*0.6,nColours*0.75,nColours*0.9,nColours];
+nColourOut=1:nColours; % Gives a nice continuous palette.
+%                          DR    R     A    Y    G     B
+cRed=interp1(nColourIn,  [0.62, 0.9, 1    , 1  , 0  , 0.46],nColourOut);
+cGreen=interp1(nColourIn,[0   , 0  , 0.52 , 1  , 0.8, 0.63],nColourOut);
+cBlue=interp1(nColourIn, [0.2 , 0.2, 0    , 0  , 0  , 0.83],nColourOut);
+colourSpec=[cRed;cGreen;cBlue]';
+
+if isfield(plotOPTS,'altColours') && plotOPTS.altColours==1
+    clear nColours nColourIn nColourOut colourSpec
+    % Build a colour palette which matches Jerry's ranges.
+    % Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
+    nColours=200;
+    nColourIn=[1,nColours*0.1,nColours*0.2,nColours*0.3,nColours*0.4,nColours];
+    nColourOut=1:nColours; % Gives a nice continuous palette.
+    %                          DB    B     LB   DG    G     LG
+    cRed=interp1(nColourIn,  [0   , 0   , 0   , 0.04, 0 , 0.2 ],nColourOut);
+    cGreen=interp1(nColourIn,[0   , 0   , 0.2 , 0.51, 0.6, 0.76 ],nColourOut);
+    cBlue=interp1(nColourIn, [0.4 , 0.6 , 0.79, 0.78, 0.6  , 0],nColourOut);
+    colourSpec=flipud([cRed;cGreen;cBlue]');
+end
+
+dataToUse=single(FVCOM.(plotOPTS.var_plot));
+
+% 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
+% Check if we're doing the cumulative difference. For this to work, you
+% need to start your data import before the start of the leak.
+% if isfield(plotOPTS,'summed_ph') && plotOPTS.summed_ph
+%     phDiff=cumsum(diff(phData,[],2));
+% else
+%     % Otherwise, just do the difference from the start, ignoring all
+%     % previous steps i.e. for each successive time step, calculate the
+%     % difference between the current time step and the background
+%     % level.
+%     phDiff=phData-repmat(bgPH,1,size(phData,2));
+% end
+
+% Check if we're doing the cumulative difference
+% if isfield(plotOPTS,'summed_ph') && plotOPTS.summed_ph
+%     phDiff=cumsum(phDiff,2);
+% end
+
+% Since we're not doing the cumulative difference, just do the difference
+% relative to the background value
+phDiff=phData-repmat(bgPH,1,size(phData,2));