How to run a function through multiple arrays

I have a function that requires input from arrays representing data from around the globe (see below). The function is designed to solve for one value at a time, but I need to solve for all the points in my global data set (i.e. I have arrays for temp, s, alk, dic that correspond to specific locations around the globe). What is the best way to loop this function through all the arrays?

function[pco2,pH,co2,hco3,co3] = co3eq(temp,s,z,alk,dic)
% Function to calculate pCO2 (microatm), pH, H2CO3* (micromol/kg), bicarbonate ion (micromol/kg), and
% carbonate ion (micromol/kg) concentrations from Alkalinity (micromol/kg) and
% DIC (micromol/kg).  This function requires temperature (degrees C),
% salinity (ppt), and depth (m), and uses the equations in Zeebe and
% Wolf-Gladrow (2000).

t = temp + 273.15;
Pr = z/10;
alk = alk*.000001;
dic = dic*.000001;
R = 83.131;

% Calculate total borate from chlorinity
tbor = .000416 * s / 35.0;

% Calculate Henry's Law coeffieicent, K0 (Weiss, 1974)
U1 = -60.2409 + 93.4517 * (100/t) + 23.3585*log(t/100); % note that *log* in matlab is the natural log, ln
U2 = s * (.023517 - .023656 * (t/100) + .0047036 * (t/100)^2);
KH = exp(U1+U2);

% Calculate KB from temp and salinity (Dickson, 1990)
KB = exp((-8966.9 - 2890.53 * s^0.5 - 77.942 * s + 1.728 * s^1.5 - 0.0996 * s^2)/t ...
    + 148.0248 + 137.1942 * s^0.5 + 1.62142 * s - (24.4344 + 25.085 * s^0.5 + ...
    0.2474 * s) * log(t) + 0.053105 * s^0.5 * t);

% Calculate K1 and K2 (Luecker et al., 2000)
K1 = 10^(-(3633.86/t - 61.2172 + 9.67770 * log(t) - 0.011555 * s + 0.0001152 * s^2));
K2 = 10^(-(471.78/t + 25.92990 - 3.16967 * log(t) - 0.01781*s + 0.0001122 * s^2));

% Pressure variation of K1, K2, and KB (Millero, 1995)
dvB = -29.48 + 0.1622 * temp - .002608 * (temp)^2;
dv1 = -25.50 + 0.1271 * temp;
dv2 = -15.82 - 0.0219 * temp;
dkB = -.00284;
dk1 = -.00308 + 0.0000877 * temp;
dk2 = .00113 - .0001475 * temp;
KB = (exp(-(dvB / (R * t)) * Pr + (0.5 * dkB / (R * t)) * Pr^2)) * KB;
K1 = (exp(-(dv1 / (R * t)) * Pr + (0.5 * dk1 / (R * t)) * Pr^2)) * K1;
K2 = (exp(-(dv2 / (R * t)) * Pr + (0.5 * dk2 / (R * t)) * Pr^2)) * K2;

% temperature dependence of KW (DoE, 1994)
KW1 = 148.96502 - 13847.26 / t - 23.65218 * log(t);
KW2 = (118.67 / t - 5.977 + 1.0495 * log(t)) * s^.5 - 0.01615 * s;
KW = exp(KW1 + KW2);

% solve for H ion (Zeebe and Wolf-Gladrow, 2000)
a1 = 1;
a2 = (alk + KB + K1);
a3 = (alk * KB - KB * tbor - KW + alk * K1 + K1 * KB + K1 * K2 - dic * K1);
a4 = (-KW  * KB + alk * KB * K1 - KB * tbor * K1 - KW * K1 + alk * K1 * K2 ...
    + KB * K1 * K2 - dic * KB * K1 - 2 * dic * K1 * K2);
a5 = (-KW * KB * K1 + alk * KB * K1 * K2 - KW * K1 * K2 - KB * tbor * K1 * K2 ...
    - 2 * dic * KB * K1 * K2);
a6 = -KB * KW * K1 * K2;
p = [a1 a2 a3 a4 a5 a6];
r = roots(p);
h = max(real(r));

% calculate bicarbonate, carbonate, and aqueous CO2 using DIC, Alk, and H+
format short g;
hco3 = dic/ (1 + h/K1 + K2/h) * 1e6;
co3 = dic / (1 + h/K2 + h * h / (K1 * K2)) * 1e6;
co2 = dic / (1 + K1/h + K1 * K2 / (h * h)) * 1e6;
pco2 = co2 / KH ;
pH = -log10(h);

% calculate B(OH)4 and OH
BOH4 = KB * tbor / (h + KB);
OH = KW / h;

% recalculate DIC and Alk to check calculations
Ct = (hco3 + co3 + co2) * 1e6;
At = (hco3 + 2*co3 + BOH4 + OH - h) * 1e6;

Answer 1

[Walter Roberson]

Everything before the roots() can be vectorized, slightly rewritten to be able to work on arrays of values (such as might be produced with ndgrid). However, to find the roots of an array of numbers, you are going to need to use arrayfun() with 'Uniform', 0. Once you have that, cellfun to find the max of the real parts of the roots, and you are back to a grid of values, and can proceed in vectorized manner.

[Brigitta Rongstad]

Just tried this (thank you by the way--silly that I didn't try this before). However, I'm still getting an error from the line of the roots function. Should I be setting up the arrayfun code differently?

Error using arrayfun
All of the input arguments must be of the same size and shape.
Previous inputs had size 180 in dimension 1. Input #3 has size 1

Error in co3eq (line 59)
[r] = arrayfun(@roots,p,'UniformOutput','false');

[Walter Roberson]

[r] = arrayfun(@roots,p,'UniformOutput',false);

However, that is not really what you want. You want

r = arrayfun(@A1,A2,A3,A4,A5,A6) roots([A1,A2,A3,A4,A5,A6]), a1, a2, a3, a4, a5, a6, 'Uniform', 0);

Your arrays are probably 2D, so when you used

p = [a1, a2, a3, a4, a5, a6]

then you got out a wide 2D array, and then your arrayfun was iterating over each element of that 2D array. If you had used

p = cat(3, a1, a2, a3, a4, a5, a6);

then you could at least worked a stack at a time, roots(squeeze(p(J,K,:)) but that is not that convenient. The work-arounds to be able to nicely do one "element" at a time such as by way of cellfun are more costly then the above code.