% This matlab function calls the fortran subroutine bulk_v2d.f 
% The matlab executable version of bulk_v2d.f (bulk_v2d.mex) is created 
% using the gateway routine bulk_v2dg.f. See the header comments in bulk_v2dg.f
% for instruction on how to compile and envoke bulk_v2d.mex
%
% The advantage of using bulk_v2d_drv.m  over the mex script is that you can 
% get help, and you can use optional inputs and outputs.
%                                
%                               
% bulk_v2d.f is adapted from bulk flux version 2.0 
%....................................................................
%  
%  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)
%     jwarm (=2 in first line; =0 if SST by IR radiometer)
%     time (time - fraction of a day)
%
%  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 m**2/s**2  
%     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
%     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 delta T parameter
%     dter  - cool skin parameter
%......................................................................
%    
% [QH,QE,RF,TAU,ws_h,Ta_h,qq_h,T0,tstar,...
%       qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter]...
%    =bulk_v2d_drv_matlab(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,rs,...
%                                     rl,rain,jcool,jwarm,time);

  function [QH,QE,RF,TAU,ws_h,Ta_h,qq_h,T0,tstar,...
          qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter]=...
bulk_v2d_drv_matlab(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,...
                           rs,rl,rain,jcool,jwarm,time);
             
if (nargin == 10);
    zi = 600*ones(size(sst));
    jcool=zeros(size(sst));
    jwarm=zeros(size(sst));
    time=zeros(size(sst));
    rain=zeros(length(sst));
    rs=zeros(length(sst));
    rl=zeros(length(sst));
'ok'
end
if(length(zi)==1)
zi=zi*ones(size(sst));
end
if(jwarm(1)==2)
if(length(jwarm)==1)
jwarm=ones(size(sst));
jwarm(1)=2;
end
else
jwarm=zeros(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 (jcool(1)==1)
if(length(jcool)==1)
jcool=ones(size(sst));
end
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
[QH,QE,RF,TAU,ws_h,Ta_h,qq_h,T0,...
tstar,qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter]...
=bulk_v2d_matlab(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,zi,...
            rs,rl,rain,jcool,jwarm,time);
          
jcool=jcool(1);jwarm=jwarm(1);hUm=hUm(1);hTm=hTm(1);hTs=hTs(1);
ts_depth=ts_depth(1);