function [varargout] = gregg90_(DATE,LATITUDE,LONGITUDE,P,Hoz,WV,AM,WM,W,RH,V)

% GREGG90_ - Simple spectral solar irradiance model 
%           for cloudless maritime atmospheres.
%
% Use As: gregg90_(DATE,LATITUDE,LONGITUDE,P,Hoz,WV,AM,WM,W,RH,V,writeopt)
% Inputs: DATE      = GMT Date (YYMMDD.HHMMSS)
%         LATITUDE  = Latitude (N+/S-  DD.MMmm)
%         LONGITUDE = Longitude (W-/E+ DD.MMmm)
%         P         = Atmospheric Pressure (mb)
%         Hoz       = Ozone Scale Hieght (cm)
%         WV        = Total Precipitable Water (cm)
%         AM        = Air Mass Type (1-10)
%         WM        = Windspeed averaged over previous
%                     24 hours (m/s)
%         W         = Instantaneous Wind speed (m/s)
%         RH        = Relative Humidity (%)
%         V         = Visual Range (km)'
%         writeopt  = Save the file to disk (0 = No (default), 1 = Yes)
%         AEROmd    = Which aerosol model to use
%                     1 = Gregg and Carder's (default)
%                     2 = Gathmans
% Output: varargout = Irradiances (uW/cm^2/nm) and Transmittances
%
% Example: gregg90_(980310.2241,21,157,1013,[],2,1,2,1,80,20)
%
% Ref: Gregg, W.W. and K.L. Carder. 1990. A Simple spectral
%      solar irradiance model for cloudless maritime atmospheres.
%      Limnology and Oceanography. Vol. 35(8): 1657-1675.

% 1 Jun 99; S. Flora - translated from Gregg and Carder, 1990

if nargin == 0
   help gregg90_
   return
end   

if nargin < 11
  disp('  GREGG90_ Error (11 Inputs Required)')
  return
end   

if nargin < 5
  Hoz = [];
end  
 
[Year,Month,Day,Hour,Minute,Second,JDAY,GMT] = timeupk_(DATE,0);  % Unpack the DATE
if Year < 30
  Year = Year + 2000;
else
  Year = Year + 1900;
end

%.............................................................................................
% Read in the File Containing Gregg and Carters Data.
mlglobal
global prntopt
prntopt = 0;
load('gregg90.mat');	            									% Read in the wavelength Ho, aoz, aw and ao
lambda = Xdat(:,1);				            									% Wavelengths (330-700 nm)
Ho 	   = Xdat(:,2);     	              					      % Mean Extraterestrial irradiance (uW/cm^2/nm)
aoz    = Xdat(:,3);				            									% Ozone absorbtion coefficients (1/cm)
aw     = Xdat(:,4);				            									% Water absorbtion coefficients (1/cm)
ao     = Xdat(:,5);			                      					% Oxygen absorbtion coefficients (1/cm)

%.............................................................................................
% Extraterrestrial solar irradiance (Eq 12, pg 1661)

e = 0.0167;											                        % Earth's orbital eccentricity
Fo = Ho.*(1+e*cos(2.*pi*(JDAY-3)/365))^2;             	% Extraterrestrial irradiance corrected for earth-sun distance

%.............................................................................................
% Calculate the zenith angle from Date.Time and Latitude and longitude

[DECL,GHA,DIST,EQ_TIME] = almanac_(Year,Month,Day,GMT);  % Moderate precision solar ephemeris
LATDD = dms2deg_(LATITUDE,0);                            % Convert to DD.dd
LONDD = dms2deg_(LONGITUDE,0);                           % Convert to DD.dd
[ALT,AZ] = altazm_(DECL,GHA,LATDD,-LONDD);                % Calculate azimuth and altitude
ZEN = (90 - ALT);                                        % Zenith angle in degrees
ZENR = (2*pi/360)*(ZEN);                                 % Zenith angle in radians
%.............................................................................................
% Atmospheric Path Length (Eq 13 & 14, pg 1661)

M   = 1/(cos(ZENR) + 0.15*((93.885 - ZEN)^(-1.253)));    % Atmosphere Path Length (Eq. 13)

Moz = 1.0035/(((cos(ZENR).^2) + 0.007).^(.5));           % Ozone Atmosphere Path Length (Eq. 14)

%.............................................................................................
% Rayleigh scattering (Eq 15 and 16, pg 1661)

Po = 1013.25;                                            % Standard Atmospheric Pressure (mb)
Mprime = M.*(P./Po);                                     % Path Length Corrected for nonstandard pressure

% CAUTION lambda are in micrometers not nanmeters
Tr = exp(-Mprime./(115.6406.*((lambda./1000).^4) - 1.335.*((lambda./1000).^2))); % Rayleigh Scattering

%.............................................................................................
% Ozone Absorbtion (Eq 17, pg 1661)

Toz = exp(-aoz.*Hoz.*Moz);                               % Ozone Absorbtion

%.............................................................................................
% Oxygen and Water Vapor Absorbtion (Eq 18 and 19, pg 1661)

To = exp((-1.4.*ao.*Mprime)./((1 + 118.3.*ao.*Mprime).^0.45));   % Oxygen Absorbtion

Tw = exp((-0.2385.*aw.*WV.*M)./((1 + 20.07.*aw.*WV.*M).^0.45));  % Water Absorbtion

%.............................................................................................
% Aerosol scattering and absorbtion (Eq 20-36, pg 1662-1666)

% Gregg and Carder version of A1-A3
R = 0.05;                                                % Factor to correct for overestimating the large particle amplitude
A(1) = 2000.*(AM.^2);                                    % Component 1 of amplitude function (eq 21), continental-derived component
A(2) = max([5.866.*(WM - 2.2) 0.5]);                     % Component 2 of amplitude function (eq 22), equilibrium sea-spray
A(3) = max([0.000014 R.*0.01527.*(W - 2.2)]);               % Component 3 of amplitude function (eq 23), ephemeral component

%      1st   2nd    3rd  Component
ro  = [0.03  0.24   2.0];                               % Mode radius of each component (um)
r   = [0.1   1.0    10 ];                               % Aerosol Particle Radius (um)

f = (((2 - (RH./100))./(6.*(1 - RH./100))).^(1/3));     % Relates particle size to relative humidity (unitless)

% Loop through the three particle components to get dN/dr for the three particle radius's
% dNdr will be used to get C and gamma.
dNdr(1:3) = 0;
for k = 1:3
  for i = 1:3
    dNdr(k) = dNdr(k) + A(i).*exp(-((log(r(k)./(f.*ro(i)))).^2)./f);
  end
end

[p, s] = polyfit(log10(r),log10(dNdr),1);
gamma = p(1);                                           % Slope = gamma (The Junge exponent)
C = p(2);                                               % Offset = C    (Amplitude function)

alpha = (-1)*(gamma + 3);                               % Angstrom exponent 

ca550 = 3.91/V;                                         % Aerosol total extinction at 550 nm
Ha    = 1;                                              % Aerosol Scale Height (km) - assumed
Ta550 = ca550*Ha;                                       % Aerosol Optical thickness at 550 nm
beta  = Ta550./((550/1000).^((-1)*alpha));              % Turbidity Coefficient representing the aerosol consentration

taua  = beta.*((lambda./1000).^((-1)*alpha));           % Aerosol Optical Thickness 

Ta    = exp((-1).*taua.*M);                             % Transmitance after aerosol scattering

% Two other aersol characteristics are required, Fa and wa for radiative transfer calculations
cos0_ = -0.1417.*alpha + 0.82;                          % Asymmetry parameter (a function of the aerosol size distribution)
if alpha < 0
  cos0_ = 0.82;                               
elseif alpha > 1.2
  cos0_ = 0.65;
end

B3 = log(1 - cos0_);
B2 = B3.*(0.0783 + B3.*(-0.3824 - 0.5874.*B3));
B1 = B3.*(1.459 + B3.*(0.1595 + 0.4129.*B3));
Fa = 1 - 0.5.*exp(cos(ZENR).*(B1 + B2.*cos(ZENR)));    % Anisotropy factor for the aerosol scattering phase function

wa = (-0.0032.*AM + 0.972).*exp((3.06E-4).*RH);        % Single scattering albedo
%.............................................................................................
% Surface Reflectance (Eq 37-47, pg 1666-1667)

pa = 1.2E3;                                            % Density of air (g/m^3)
D1 = 2.2E-5;                                           % Coeffecient 1 relating windstress to foam reflectance
D2 = 4.0E-4;                                           % Coeffecient 2 relating windstress to foam reflectance
D3 = 4.5E-5;                                           % Coeffecient 3 relating windstress to foam reflectance
D4 = 4.0E-5;                                           % Coeffecient 4 relating windstress to foam reflectance
if W <= 4
  CD = NaN;
  pf = 0;                                              % Surface reflectance due to Sea Foam
elseif W <= 7
  CD = (0.62 + 1.56.*W.^(-1))*10^-3;                   % Drag Coefficient 
  pf = D1.*pa.*CD.*W.^2 - D2;                          % Surface reflectance due to Sea Foam
else 
  CD = (0.49 + 0.065.*W)*10^-3;                        % Drag Coefficient 
  pf = (D3.*pa.*CD - D4).*W.^2;                        % Surface reflectance due to Sea Foam
end

ZENr = NaN;
if W <= 2 | abs(ZEN) <= 40                          
  nw = 1.341;                                          % Index of refraction of seawater
  ZENr = asin(sin(ZENR)/nw);                           % The refracted angle (radians) ??????????
  pd1  = ((sin(ZENR-ZENr)^2)./(sin(ZENR+ZENr)^2));
  pd2  = ((tan(ZENR-ZENr)^2)./(tan(ZENR+ZENr)^2));
  pd   = 0.5*(.5*pd1 +pd2);                            % (Fresnel) Specular reflectance for direct irradiance (~flat sea or high zenith)
else
  b    = (-7.14e10-4)*W + 0.0618;
  pdsp = 0.0253*exp(b*(ZEN - 40));                     % Direct Specular reflectance for direct irradiance (rough sea and low zenith)
  pd   = pdsp + pf;                                    % Specular reflectance for direct irradiance 
end

if W <= 4
  pssp = 0.066;                                        % Diffuse Specular reflectance for direct irradiance
elseif W > 4
  pssp = 0.057;                                        % Diffuse Specular reflectance for direct irradiance
end
ps = pssp + pf;                                        % Specular reflectance for diffuse irradiance
%.............................................................................................
% The Model (Eq 9-11, pg 1660) - PUTTING IT ALL TOGETHER
Edd = Fo.*cos(ZENR).*Tr.*Ta.*Toz.*To.*Tw;                   % Direct Irradiance 
Taa = exp((-1).*(1 - wa).*M.*taua);                         % Transmittance after aerosol absorbtion
Ir  = Fo.*cos(ZENR).*Toz.*To.*Tw.*Taa.*(1-Tr.^0.95).*(0.5); % Diffuse Component arising from Rayleigh Scattering
Tas = exp((-1).*wa.*taua.*M);                               % Transmittance due to aerosol scattering
Ia  = Fo.*cos(ZENR).*Toz.*To.*Tw.*Taa.*Tr.^1.5.*(1-Tas)*Fa; % Diffuse Component arising from Aerosol Scattering
Eds = (Ir + Ia);                                            % Diffuse Irradiance 

Ed = Edd.*(1 - pd) + Eds.*(1 - ps);                         % Downwelled Irradiance Just below the Surface
