function [ q900, q600 ] = qs900b( T600, vmax )
% Calculates saturation specific humidity at 600 hPa and saturation
% specific humidity at pressure 'pref' given 600 hPa T (K) and assuming a
% moist adiabatic lapse rate.  
%
%  This version is faster than qs900 and takes into account increased
%  eyewall entropy, which requires input of circular wind. 
%
% Copyright WindRiskTech, LLC  October, 2017
%-----------------------------------------------------------------------
pref=950;    %  Pressure to find qs at
%
cp=1005;
Rv=491;
Rd=287;
Lv=2.5e6;
%
c1=Lv./Rv;
c2=Rd*log(pref/600);
c3=1.6/100;  %  The factor 1.6 rescales the square of the maximum wind speed 
             %  to account for surface reduction from gradient wind. 100=T_s-T_t. 
vmax=vmax.*1852/3600; % Convert from knots to m/s.
Tc=T600-273.15;
Tc=max(Tc,-50);
es=6.112.*exp(17.67.*Tc./(243.5+Tc));
q600=0.622.*es./(600-es);
%
clear es Tc
T=T600+20;   % First guess
[qn,qm]=size(q600);
q900=zeros(qn,qm);
%
for i=1:qn
    for j=1:qm
        if T600(i,j) > 100 % Exclude zero fill values of T600
             for it=1:5
                 Tc=T(i,j)-273.15;
                 es=6.122.*exp(17.67.*Tc./(243.5+Tc));
                 qs=0.622*es./(pref-es);
                 er=cp*log(T(i,j)./T600(i,j))+Lv.*(qs./T(i,j)-q600(i,j)./T600(i,j))-c2-c3.*vmax(i,j).^2;
                 derdT=(cp*T(i,j)+Lv.*qs.*(c1./T(i,j)-1))./T(i,j).^2;
                 T(i,j)=T(i,j)-er./derdT;
             end
             q900(i,j)=qs;
        end
    end
end
%
end