function [mystruct_out]=gmaze_bulk_coare(mystruct_in,varargin);

%varargin       [hu ht hq]
%               [hu ht hq],albedo
%       [albedo HAS AN impact here]

%version with shortened iteration; modified Rt and Rq
%uses wave information wave period in s and wave ht in m
%no wave, standard coare 2.6 charnock:  jwave=0 
%Oost et al.  zo=50/2/pi L (u*/c)^4.5 if jwave=1
%taylor and yelland  zo=1200 h*(L/h)^4.5 jwave=2

%inputs:
atemp_in=mystruct_in.atemp;
aqh_in=mystruct_in.aqh;
uwind_in=mystruct_in.uwind;
vwind_in=mystruct_in.vwind;
swdn_in=mystruct_in.swdn;
lwdn_in=mystruct_in.lwdn;
slp_in=mystruct_in.slp;
precip_in=mystruct_in.precip;
sst_in=mystruct_in.sst;

%compute wind speed:
wspeed=uwind_in.*uwind_in+vwind_in.*vwind_in;
wspeed=sqrt(wspeed); %wspeed=max(wspeed,0.5);

%compute bulk water spec hum:
nlon=size(wspeed,1); nlat=size(wspeed,2);
Qs=reshape(qsee([reshape(sst_in,[nlon*nlat 1]) reshape(slp_in,[nlon*nlat 1])]),[nlon nlat]);
%	=> gives qsat in g/kg
%in ... exf_bulk_largeyeager04:
%cvapor_fac=640380;
%cvapor_exp=5107.4;
%saltsat=0.980;
%atmrho=1.2;
%Qs = cvapor_fac*exp(-cvapor_exp./(sst_in+273.16))*saltsat/atmrho; 
%	=> gives qsat in kg/kg


%change name of variable:
u      = wspeed;        	% wind speed (m/s)  at height zu (m)
us     = 0*ones(size(u));      	% surface current speed in the wind direction (m/s)
ts     = sst_in;		% bulk water temperature (C) if jcool   = 1, interface water T if jcool   = 0  
t      = atemp_in-273.16; 	% bulk air temperature (C), height zt
Qs     = Qs/1000;%bulk water spec hum out of qsee -> /1000 => kg/kg
Q      = aqh_in; %bulk air spec hum is in kg/kg (from input)
Rs     = swdn_in;		% downward solar flux (W/m^2)
Rl     = lwdn_in;      		% downard IR flux (W/m^2)
rain   = precip_in/3600;	% rain rate (mm/hr)
zi     = 600*ones(size(u));	% PBL depth (m)
P      = slp_in;		% Atmos surface pressure (mb)
grav   = 9.81*ones(size(u));	% grv(lat)
if nargin==4
albed  = varargin{2}*ones(size(u));     % albed
else
albed  = 0.1*ones(size(u));	% albed
end

% Vectorization:
[ny nx] = size(u);
u  = u(:);
us = us(:);
ts = ts(:);
t  = t(:);
Qs = Qs(:);
Q  = Q(:);
Rs = Rs(:);
Rl = Rl(:);
rain = rain(:);
zi = zi(:);
P  = P(:);
grav = grav(:);
albed = albed(:);

% 1D fields:
if nargin==1
zu     = 10;  % wind speed measurement height (m)
zt     = 10;  % air T measurement height (m)
zq     = 10;  % air q measurement height (m) (specific humidity in g/kg)
else
zu=varargin{1}(1);
zt=varargin{1}(2);
zq=varargin{1}(3);
end
jcool  = 1; % implement cool calculation skin switch, 0 = no, 1 = yes
jwave  = 0; % implement wave dependent roughness model
twave  = 0; % wave period (s)
hwave  = 0; % wave height (m)
% add the adjusted height
zsu    = 10; %wind speed output height (m)
zst    = 10; %air T output height (m)
zsq    = 10; %air q output height (m)


% BEGIN THE COMPUTATION:


     %******.*.*.*.*.*   set constants .*.*.*.*.*.*.*.*.*.*.*.*.*
     Beta = 1.2;
     von  = 0.4;
     fdg  = 1.00;
     tdk  = 273.16;
     Rgas = 287.1;
     cpa  = 1004.67;
     be   = 0.026;
     cpw  = 4000;
     rhow = 1022;
     visw = 1e-6;
     tcw  = 0.6;
     dter = 0.3;

     %.*.*.*.*.*.*.*.*.*.*.*.*.*  air constants .*.*.*.*.*.*.*.*.*.*.*.*
     Rgas = 287.1;
     Le   = (2.501-.00237.*ts).*1e6;
     cpa  = 1004.67;
     cpv  = cpa.*(1+0.84.*Q);
     rhoa = P.*100./(Rgas.*(t+tdk).*(1+0.61.*Q));
     visa = 1.326e-5.*(1+6.542e-3.*t+8.301e-6.*t.*t-4.84e-9.*t.*t.*t);

     %.*.*.*.*.*.*.*.*.*.*.*.*  cool skin constants  .*.*.*.*.*.*.*
     Al   = 2.1e-5.*(ts+3.2).^0.79;
     be   = 0.026;
     cpw  = 4000;
     rhow = 1022;
     visw = 1e-6;
     tcw  = 0.6;
     bigc = 16.*grav.*cpw.*(rhow.*visw).^3./(tcw.*tcw.*rhoa.*rhoa);
     wetc = 0.622.*Le.*Qs./(Rgas.*(ts+tdk).^2);

     %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*   wave parameters  .*.*.*.*.*.*.*.*.*
     lwave=grav./2./pi.*twave.^2;
     cwave=grav./2./pi.*twave;

     %.*.*.*.*.*.*.*.*.*.*.*.*.*.*  compute aux stuff .*.*.*.*.*.*.*
     % Net short wave
     Rns = Rs.*(1-albed);
     % Net long wave   
     Rnl = 0.97.*(5.67e-8.*(ts-0.3.*jcool+tdk).^4-Rl);

     %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*   Begin bulk loop .*.*.*.*.*.*.*

     %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*  first guess .*.*.*.*.*.*.*.*.*.*.*.*
     du=u-us;
     dt=ts-t-.0098.*zt;
     dq=Qs-Q;
     ta=t+tdk;
     ug=.5;
     dter=0.3; 
     dqer=wetc.*dter;
     ut=sqrt(du.*du+ug.*ug);
     % whoi using zsu (could be 10 or 2m)instad of 10m
     %u10=ut.*log(10./1e-4)./log(zu./1e-4);
     u10=ut.*log(zsu./1e-4)./log(zu./1e-4);
     usr=.035.*u10;
     zo10=0.011.*usr.*usr./grav+0.11.*visa./usr;
     Cd10=(von./log(zsu./zo10)).^2;
     Ch10=0.00115;
     Ct10=Ch10./sqrt(Cd10);
     zot10=zsu./exp(von./Ct10);
     Cd=(von./log(zu./zo10)).^2;
     Ct=von./log(zt./zot10);
     CC=von.*Ct./Cd;
     Ribcu=-zu./zi./.004./Beta.^3;
     Ribu=-grav.*zu./ta.*((dt-dter.*jcool)+.61.*ta.*dq)./ut.^2;
     nits=3;
     if Ribu<0;
	zetu=CC.*Ribu./(1+Ribu./Ribcu);
     else;
	zetu=CC.*Ribu.*(1+27./9.*Ribu./CC);
     end;		
     L10=zu./zetu;
     if zetu>50;
	nits=1;
     end;
     usr=ut.*von./(log(zu./zo10)-psiu_30(zu./L10));
     tsr=-(dt-dter.*jcool).*von.*fdg./(log(zt./zot10)-psit_30(zt./L10));
     qsr=-(dq-wetc.*dter.*jcool).*von.*fdg./(log(zq./zot10)-psit_30(zq./L10));

     tkt=.001;
	
     charn=0.011;
     if ut>10
        charn=0.011+(ut-10)./(18-10).*(0.018-0.011);
     end;
     if ut>18
        charn=0.018;
     end;

     %disp(usr)

     %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*  bulk loop .*.*.*.*.*.*.*.*.*.*.*.*
     for i=1:nits;

        zet=von.*grav.*zu./ta.*(tsr.*(1+0.61.*Q)+.61.*ta.*qsr)./(usr.*usr)./(1+0.61.*Q);
        if jwave==0;zo=charn.*usr.*usr./grav+0.11.*visa./usr;end;
        if jwave==1;zo=50./2./pi.*lwave.*(usr./cwave).^4.5+0.11.*visa./usr;end;%Oost et al
        if jwave==2;zo=1200.*hwave.*(hwave./lwave).^4.5+0.11.*visa./usr;end;%Taylor and Yelland
        rr=zo.*usr./visa; % roughness Reynolds number
        L=zu./zet;
        zoq=min(1.15e-4,5.5e-5./rr.^.6);
        zot=zoq;
        usr=ut.*von./(log(zu./zo)-psiu_30(zu./L));
        tsr=-(dt-dter.*jcool).*von.*fdg./(log(zt./zot)-psit_30(zt./L));
        qsr=-(dq-wetc.*dter.*jcool).*von.*fdg./(log(zq./zoq)-psit_30(zq./L));
        Bf=-grav./ta.*usr.*(tsr+.61.*ta.*qsr);
        if Bf>0
           ug=Beta.*(Bf.*zi).^.333;
        else
           ug=.2;
        end;
        ut=sqrt(du.*du+ug.*ug);
        Rnl=0.97.*(5.67e-8.*(ts-dter.*jcool+tdk).^4-Rl);
        hsb=-rhoa.*cpa.*usr.*tsr;
        hlb=-rhoa.*Le.*usr.*qsr;
        qout=Rnl+hsb+hlb;
        dels=Rns.*(.065+11.*tkt-6.6e-5./tkt.*(1-exp(-tkt./8.0e-4))); 	% Eq.16 Shortwave
        qcol=qout-dels;
        alq=Al.*qcol+be.*hlb.*cpw./Le;					% Eq. 7 Buoy flux water

        if alq>0;
           xlamx=6./(1+(bigc.*alq./usr.^4).^.75).^.333;			% Eq 13 Saunders
           tkt=xlamx.*visw./(sqrt(rhoa./rhow).*usr);			%Eq.11 Sub. thk
        else
           xlamx=6.0;
           tkt=min(.01,xlamx.*visw./(sqrt(rhoa./rhow).*usr));		%Eq.11 Sub. thk
        end;

        dter=qcol.*tkt./tcw;%  Eq.12 Cool skin
        dqer=wetc.*dter;

    end;%bulk iter loop

    tau =  rhoa.*usr.*usr.*du./ut;                %stress
    hsb = -rhoa.*cpa.*usr.*tsr;
    hlb = -rhoa.*Le.*usr.*qsr;
    Rnl = 0.97.*(5.67e-8.*(ts-dter.*jcool+tdk).^4-Rl);

    %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*   rain heat flux .*.*.*.*.*.*.*.*
    if 0     
    dwat=2.11e-5.*((t+tdk)./tdk).^1.94;%! water vapour diffusivity
    dtmp=(1.+3.309e-3.*t-1.44e-6.*t.*t).*0.02411./(rhoa.*cpa); 	%!heat diffusivity
    alfac= 1./(1+(wetc.*Le.*dwat)./(cpa.*dtmp));      	%! wet bulb factor
    RF= rain.*alfac.*cpw.*((ts-t-dter.*jcool)+(Qs-Q-dqer.*jcool).*Le./cpa)./3600;
    %.*.*.*.*.*.* add the momentum flux due to rainfall by kelan .*.*.*.*.*.*

    tau_r =0.85 .*rain./3600.*du; %%du==wind spd relative to sea surface

    %.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*   Webb et al. correection  .*.*.*.*.*.*.*.*.*.*.*.*
    wbar=1.61.*hlb./Le./(1+1.61.*Q)./rhoa+hsb./rhoa./cpa./ta;%formulation in hlb already includes webb
    %wbar=1.61.*hlb./Le./rhoa+(1+1.61.*Q).*hsb./rhoa./cpa./ta;
    hl_webb=rhoa.*wbar.*Q.*Le;
    %.*.*.*.*.*.*.*.*.*.*.*.*.*.*   compute transfer coeffs relative to ut @meas. ht .*.*.
    Cd=tau./rhoa./ut./max(.1,du);
    Ch=-usr.*tsr./ut./(dt-dter.*jcool);
    Ce=-usr.*qsr./(dq-dqer.*jcool)./ut;
    %.*.*.*.*.*.*.*.*.*.*.*.*  10-m neutral coeff realtive to ut .*.*.*.*.*.*.*.*
    Cdn_10=von.*von./log(zsu./zo)./log(zsu./zo); %von=0.4; == Frank's CD
    Chn_10=von.*von.*fdg./log(zsu./zo)./log(zsu./zot);
    Cen_10=von.*von.*fdg./log(zsu./zo)./log(zsu./zoq);
    end % 0/1    
    Cdn_10=von.*von./log(zsu./zo)./log(zsu./zo); %von=0.4; == Frank's CD


    %.*.*.*Kelan adds azeta() to adjust met variables to new heights 
    %.*.*.*copy from cor3_0af.for   
    % zsu % adjusted wind speed height
    % zst % adjusted air temp height
    % zsq % adjusted specific humidity height
    % ts - dter.*jcool === sst of Bradley's;

    if 0
    % kelan add this function call azeta.m in matlab version
    % Frank missed this call in his cor30a.for version
    zetU=azeta(von,grav,ta,Q,usr,tsr,qsr,zsu);  
    zetT=azeta(von,grav,ta,Q,usr,tsr,qsr,zst); 
	% same as above if zst=zsu
    zetQ=azeta(von,grav,ta,Q,usr,tsr,qsr,zsq); 
	% same as above if zsq=zst=zsu

    u_zs=usr./von.*(log(zsu./zo)-psiu_30(zetU));
    t_zs=(ts-dter.*jcool)+tsr./von.*(log(zst./zot)-psit_30(zetT));
    q_zs=(Qs-dqer)+qsr./von.*(log(zsq./zoq)-psit_30(zetQ)); % !kg./kg
    qa_zs=1000..*q_zs;                              % !g./kg

    e_zs=q_zs.*P./(0.62197+0.378.*q_zs);              % !mb

    %call humidity(t_zs,p,esat_zs);                % !mb
    esat_zs = (1.0007+3.46e-6.*P).*6.1121.*exp(17.502.*t_zs./(240.97+t_zs));
    rh_zs=e_zs./esat_zs;      
    end % 0/1

    % Output:
    y(1,:,:) = -reshape(hsb,[ny nx]);
    y(2,:,:) = -reshape(hlb,[ny nx]);
    y(3,:,:) = -reshape(Rnl,[ny nx]);
    y(4,:,:) =  reshape(Rns,[ny nx]);
    y(5,:,:) =  reshape(tau,[ny nx]);
    y(6,:,:) =  reshape(Cdn_10,[ny nx]);
    y(7,:,:) =  reshape(rhoa,[ny nx]);


%   compute wind stress in two directions:
    ustress =  rhoa.*usr.*usr.*uwind_in(:)./ut;
    vstress =  rhoa.*usr.*usr.*vwind_in(:)./ut;
%   compute evaporation:
    flamb=2500000;
    evap=hlb/flamb;
%   change sign and convert from kg/m^2/s to m/s via rhoConstFresh
    rhoConstFresh=999.8;
    evap = -evap/rhoConstFresh;
%   reshape:
    hlb=-reshape(hlb,[ny nx]);
    hsb=-reshape(hsb,[ny nx]);
    ustress=reshape(ustress,[ny nx]);
    vstress=reshape(vstress,[ny nx]);
    evap=-reshape(evap,[ny nx]);

%mystruct_out=struct('hl',hlb,'hs',hsb,'evap',evap,'ustress',ustress,'vstress',vstress);
mystruct_out=struct('hl',real(hlb),'hs',real(hsb),'evap',real(evap),'ustress',real(ustress),'vstress',real(vstress));


function psi=psit_30(zet)
	x=(1-15*zet).^.5;
	psik=2*log((1+x)/2);
	x=(1-34.15*zet).^.3333;
	psic=1.5*log((1+x+x.*x)/3)-sqrt(3)*atan((1+2*x)/sqrt(3))+4*atan(1)/sqrt(3);
	f=zet.*zet./(1+zet.*zet);
   psi=(1-f).*psik+f.*psic;  

   ii=find(zet>0);
if ~isempty(ii);
	%psi=-4.7*zet;
	c=min(50,.35*zet);
   psi(ii)=-((1+2/3*zet(ii)).^1.5+.6667*(zet(ii)-14.28)./exp(c(ii))+8.525);
end;


function psi=psiu_30(zet)

	x=(1-15*zet).^.25;
	psik=2*log((1+x)/2)+log((1+x.*x)/2)-2*atan(x)+2*atan(1);
	x=(1-10.15*zet).^.3333;
	psic=1.5*log((1+x+x.*x)/3)-sqrt(3)*atan((1+2*x)/sqrt(3))+4*atan(1)/sqrt(3);
	f=zet.*zet./(1+zet.*zet);
	psi=(1-f).*psik+f.*psic;                                               
   ii=find(zet>0);
   if ~isempty(ii);

	%psi(ii)=-4.7*zet(ii);
  	%c(ii)=min(50,.35*zet(ii));
   c=min(50,.35*zet);
	psi(ii)=-((1+1.0*zet(ii)).^1.0+.667*(zet(ii)-14.28)./exp(c(ii))+8.525);
	end;

function s=qsee(y)
x=y(:,1);
p=y(:,2);
es=6.112.*exp(17.502.*x./(x+240.97))*.98.*(1.0007+3.46e-6*p);
s=es*621.97./(p-.378*es);

function y=qsat(y)
x=y(:,1);%temp
p=y(:,2);%pressure
es=6.112.*exp(17.502.*x./(x+241.0)).*(1.0007+3.46e-6*p);
y=es*622./(p-.378*es);