%dicbio_main % simple model linking phytoplankton growth & calcification for % define some global variables that can be used in other subroutines % list of global variables should match across subroutines global tc press s pco2atm0 mld P % assign values to global variables (for this example will remain % constant over integration) tc=25; % temperature (deg. C) press=0; % pressure s=36.5; % salinity pco2atm0=360; % atmospheric CO2 (uatm) mld=100; % mixed layer depth (m) % give initial conditions for variables that are loaded into a vector X0 X0(1) = 0.1; % phytoplankton (umol C/kg) X0(2) = 2000; % dissolved inorganic carbon (umol/kg) X0(3) = 2350; % total alkalinity (ueq/kg) % give parameter values P(1) = 1.0; % phytoplankton growth rate (1/day) P(2) = 30.0; % carrying capacity (umol C/kg) P(3) = 0.1; % ratio of CaCO3 to Organic Carbon production (mol/mol) % set up time vector T=[0:0.5:15]'; % use a fixed time-scale to compare cases n=length(T); % integrate initial value problem forward in time [T,X]=ode15s('dicbio',T,X0); % plot model variable results figure(1) plot(T,X(:,1),'-',T,X(:,2)-X0(2),'--',T,X(:,3)-X0(3),'-.') xlabel('Time (days)') ylabel('Phyto (-), DIC anomaly (--), and TAlk anomlay (-.)') title('Initial Value Problem; Biological DIC Drawdown') % setup some diagnostic variables on chemistry pH_vec = zeros(n,1); co3_vec = zeros(n,1); pco2_vec = zeros(n,1); % loop over time steps and compute carbonate chemistry from model DIC, Alk for i=1:n dic = X(i,2); alk = X(i,3); [pco2,co2,hco3,co3,co2eq,h]=pco2_dicalk(tc,s,press,dic,alk,pco2atm0); pH_vec(i) = -log10(h); co3_vec(i) = co3; pco2_vec(i) = pco2; end % plot model chemistry results figure(2) plotyy(T,pH_vec,T,co3_vec) xlabel('Time (days)') ylabel('pH (left axis); carbonate ion (right axis)') title('Initial Value Problem; Biological DIC Drawdown') % plot model chemistry results figure(3) plot(T,pco2_vec) xlabel('Time (days)') ylabel('pCO2') title('Initial Value Problem; Biological DIC Drawdown')