function out = ERA40atm(site_lat,site_lon,site_elv)
% Looks up mean sea level pressure and mean 1000 mb temp from ERA-40 reanalysis
% and calculates site atmospheric pressures using these as inputs to the
% standard atmosphere equation.
%
% Syntax: pressure = ERA40(site_lat,site_lon,site_elv);
%
% Requires:
% site_lat: latitude (DD). Southern hemisphere is negative.
% site_lon: longitude (DD). Western hemisphere is negative.
% Tries to deal with 0-360 longitudes gracefully.
% site_elv: elevation (m).
%
% Returns site pressure in hPa.
%
% Vectorized. Send vectors of equal length. Will fail with mystery error
% message if you don't.
%
% Note: this must load the data file ERA40.mat whenever called.
% Repeated calls to this function will be slow for this reason.
%
% Also: This function is probably OK for Antarctica, but you should
% consider using antatm.m instead.
%
% Finally: it is often better to estimate the average pressure at your
% site using a pressure-altitude relation obtained from nearby station
% data.
%
% Written by Greg Balco -- UW Cosmogenic Nuclide Lab
% balcs@u.washington.edu
% October, 2007
% Part of the CRONUS-Earth online calculators:
% http://hess.ess.washington.edu/math
%
% Copyright 2001-2007, University of Washington
% All rights reserved
% Developed in part with funding from the National Science Foundation.
%
% Modified by Nat Lifton -- Purdue University to use ERA-40 instead of
% NCEP, and also to use a variable lapse rate.
% March 2011
%
% Greg Balco checked Dec 2014.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License, version 2,
% as published by the Free Software Foundation (www.fsf.org).
% correct negative longitudes
sz = find(site_lon < 0);
site_lon(sz) = site_lon(sz) + 360;
% data load
load ERA40;
% Interpolate sea level pressure and 1000-mb temperature
% from global reanalysis data grids.
% site_T in K, site_P in hPa
site_slp = interp2(ERA40lon,ERA40lat,meanP,site_lon,site_lat);
site_T = interp2(ERA40lon,ERA40lat,meanT,site_lon,site_lat);
% site_slp = 1013.25;
% site_T = 288.15;
%
% site_T_degK = site_T + 273.15;
% More parameters
gmr = -0.03417; % Assorted constants
% dtdz = 0.0065; % Lapse rate from standard atmosphere
% Lifton Lapse Rate Fit to COSPAR CIRA-86 <10 km altitude
lr = [-6.1517E-03 -3.1831E-06 -1.5014E-07 1.8097E-09 1.1791E-10 ...
-6.5359E-14 -9.5209E-15];
dtdz = lr(1) + lr(2).*site_lat + lr(3).*site_lat.^2 ...
+ lr(4).*site_lat.^3 + lr(5).*site_lat.^4 + lr(6).* site_lat.^5 ...
+ lr(7).*site_lat.^6;
dtdz = -dtdz;
% Variable lapse rate experiments -- attempts to make lapse rate a
% physically reasonable function of temperature. No guarantees on the
% correctness of this part.
%
% Not used in final version. Probably more physically correct but
% has a very limited effect on the overall results. Code retained here as
% comments for those who are interested. If you are an atmospheric
% scientist, feel free to suggest a better way to do this.
%
% Temp assumed for lapse rate is 15 degrees less than real temp.
% Chosen because in stdatm, std T is 15 and LR is MALR for 0 deg C.
% Not sure why this is the case. Ask an atmospheric scientist.
%
%Tlr = site_T_degK - 15;
%
% Calculate the saturation vapor pressure
%
%esat = 10.*0.6112.*exp(17.67.*(Tlr-273.15)./((Tlr-273.15)+243.5));
%
% Calculate the mixing ratio and thence the lapse rate
%
%rv = 0.622.*esat./(1013.25-esat);
%Lv = 2.501e6; R = 287; E = 0.62; g = 9.8066; cpd = 1005.7;
%dtdz = g.*(1+(Lv*rv./R./Tlr))./(cpd + (Lv.^2.*rv.*E./R./(Tlr.^2)));
% Calculate site pressure using the site-specific SLP and T1000 with the
% standard atmosphere equation.
out = site_slp .* exp( (gmr./dtdz) .* ( log(site_T) - log(site_T - (site_elv.*dtdz)) ) );