Wind Stresses computation

This function computes the wind stresses from wind component

%
% USAGE: [Tx, Ty]=ra_windstr(u,v)
%
% DESCRIPTION: Function to compute wind stress from wind field data Based on Gill, 1982
% Formula and a non-linear Cd based on Large and Pond (1981), modified for low wind
% speeds (Trenberth et al., 1990)
%
% INPUTS:
% u = Zonal wind component [m/s], must be 2D
% v = Meridional wind component [m/s], must be 2D
%
% OUTPUT:
% Tx = Zonal wind stress [N/m^2]
% Ty = Meridional wind stress [N/m^2]
%
% DISCLAIMER:
% Albeit this function is designed only for academic purpose, it can be implemented in
% research. Nonetheless, author does not guarantee the accuracy.
%
% REFERENCE:
% A.E. Gill, 1982, “Atmosphere-Ocean Dynamics”, Academy Press, Vol. 30.
% W. G. Large & S. Pond., 1981,“Open Ocean Measurements in Moderate to Strong Winds”,
% J. Physical Oceanography, Vol. 11, pp. 324 - 336.
% K.E. Trenberth, W.G. Large & J.G. Olson, 1990, “The Mean Annual Cycle in Global Ocean
% Wind Stress”, J.Physical Oceanography, Vol. 20, pp. 1742 – 1760.
% ***********************************************************************************************%