% This matlab function calls the mex routine bulk_mat_v25b.mex$$ 
% The matlab executable version of bulk_mat_v25b.mex$$ is created 
% using the gateway routine bulk_mat_v25b_function.f. 
% See the header comments in bulk_v25b_function.f
% for instruction on how to compile and envoke bulk_mat_v25b.mex$$
% The exact form of the mex file is machine dependent
%
% The advantage of using bulk_v25b.m  over the mex script is that you can 
% get help, and you can use optional inputs and outputs.
%                                
%                               
% bulk_mat_v25b.f is adapted from bulk flux version 2e and 2f
% see header of bulk_mat_v25b.f for complete explanation of changes 
%....................................................................
%  
%  INPUT:
%
%     hUm (U MEASUREMENT height, m)                          scalar or vector  
%     hTm (T MEASUREMENT height, m)  assume hQm=hTm          scalar or vector
%     hUs (U STANDARD height, m)                             scalar or vector  
%     hTs (T STANDARD height, m)   assume hQs=hTs            scalar or vector 
%     ts_depth  depth of sst measurement (=0 to ignore PWP; otherwise
%                     depth in meters)
%     ws (wind speed relative to ocean speed, m/s)                vector
%     sst (sea surface temperature,  deg. C)                      vector
%     atb (air temperature, deg. C)                               vector
%     qq (specific humidity, g/kg, OR rel humidity, decimal)      vector
%     pr (pressure, mb)                                      scalar or vector
%  OPTIONAL INPUT:
%     zi (boundary-layer height, m --if N/A default is 600m)      vector 
%     rs (shortwave radiation W/m**2)                             vector
%     rl (downwelling longwave radiation W/m**2)                  vector
%     rain (average rainfall mm/hour)                             vector
%     jcool (=1 for cool skin; =0 if SST by IR radiometer)        scalar
%     jwarm (=1 for warm layer; =0 if SST by IR radiometer)       scalar
%     time (time - time of day in units of day (e.g., 0.5 == 12)
%                                             or yymmddhhmmss.ss) vector
%     lat (latitude degrees, +ve North, -ve South)
%     lon (longitude degrees, +ve East, -ve West)
%
%  OUTPUT:   (all output variables are vectors)
%
%     QH W/m**2   sensible heat flux
%     QE W/m**2   latent heat flux
%     RF W/m**2   rainfall heat flux
%     TAU N/m**2  wind stress  
%     TAUr N/m**2 stress due to rainfall
%     QE_webb W/m**2 Webb corrrection to latent heat flux
%     ws_h  - relative ws at standard height 
%     qq_h  - specific humidity or relative humidity at standard height 
%     Ta_h  - Ta at standard height 
%     T0 - skin temperature C (T0 = sst - dter + dt_wrm)
%     Ustar m/s                 
%     Qstar kg/kg
%     Tstar C
%     CD - drag coefficient
%     CH - transfer coefficient for heat
%     CE - transfer coefficient for moisture
%     RR - Roughness Reynolds number
%     RT - Roughness Reynolds number for temperature
%     RQ - Roughness Reynolds number for moisture
%     ZL - hUm/L where L is the Obukhov length
%     Z0 - roughness length
%     Zot - roughness length (temp)
%     Zoq - roughness length (humidity)
%     dt_wrm - warm layer temperature difference C
%     dter - cool skin temperature difference C
%     wg - gustiness m/s
%......................................................................
%    
% [QH,QE,RF,TAU,TAUr,QE_webb,ws_h,Ta_h,qq_h,T0,tstar,...
%       qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter,wg]...
%    =bulk_v25b(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,rs,...
%                                 rl,rain,jcool,jwarm,time,lat,lon);

  function [QH,QE,RF,TAU,TAUr,QE_webb,ws_h,Ta_h,qq_h,T0,tstar,...
          qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter,wg]=...
bulk_v25b(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,...
                           rs,rl,rain,jcool,jwarm,time,lat,lon);
             
if (nargin == 10);
    zi = 600*ones(size(sst));
    jcool=0;
    jwarm=0;
    time=zeros(size(sst),1);
    rain=zeros(length(sst),1);
    rs=zeros(length(sst),1);
    rl=zeros(length(sst),1);
    lat=zeros(length(sst),1);
    lon=zeros(length(sst),1);
end
if(length(zi)==1)
   zi=zi*ones(size(sst));
end
if (length(hUm) == 1)
   hUm = hUm*ones(size(sst));
   hTm = hTm*ones(size(sst));
end
if (length(hUs) == 1)
   hUs = hUs*ones(size(sst));
   hTs = hTs*ones(size(sst));
end
if (length(pr) == 1)
   pr = pr*ones(size(sst));
end
if (length(ts_depth) == 1)
   ts_depth = ts_depth*ones(size(sst));
end
if (length(rs)==1)
    rs=rs*ones(size(sst));
end
if (length(rl)==1)
   rl=rl*ones(size(sst));
end
if (length(rain)==1)
   rain=rain*ones(size(sst));
end
if (length(time)==1)
   time=time*ones(size(sst));
end
if (length(jcool)>1)
   jcool=jcool(1);
end
if (length(jwarm)>1)
   jwarm=jwarm(1);
end
if (length(lat)==1)
    lat=lat*ones(size(sst));
end
if (length(lon)==1)
    lon=lon*ones(size(sst));
end
%
[QH,QE,RF,TAU,TAUr,QE_webb,ws_h,Ta_h,qq_h,T0,...
tstar,qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter,wg]...
=bulk_mat_v25b(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,...
            rs,rl,rain,jcool,jwarm,time,lat,lon);