/* 'C' library for combined sea state and atmospheric boundary layer model
Developed by Mark A. Bourassa
Based on Original Bourassa-Vincent-Wood combined sea state and atmospheric
boundary layer model developed at Purdue University. Several upgrades and
options have been added while at the Center for Ocean-Atmospheric
Prediction Studies (COAPS) at Florida State University.
Dr. Bourassa is currently (Aug. 2003) an Assistant Professor at the Florida
State Universiy.
Please send comments, suggestions, and requests for added features to
Dr. Mark Bourassa
Center for Ocean-Atmospheric Prediction Studies
Florida State University
Tallahassee, FL 32306-2840
Email: bourassa@coaps.fsu.edu
Phone: (850) 644-6923
The bulk of 'in code' documentation is in the subroutine pmix_. Documentation
is also on the COAPS website: http://coaps.fsu.edu/~bourassa/bvw98_docs.html.
*/
#include "mft12.h"
/* beta-primes are either zero or unity, depending on whether or not
their respective type of roughness element exists. Beta's are the
weights for the different types of elements. */
/* global variables */
short flux_model_dum, no_capw, no_sfcten, smith88, stab_prm,
TQzo_prm, use_dh, wave_stab_intr, zo_prm;
double betac[2], betac_prime, betag[2], betag_prime, betas[2], betas_prime,
cmin, CP, cp_com, denair, hsig[2], LV, monin_inv, qstar, tstar,
u[2], u_orb[2], u_sfc[2], ustar[2], wa[2], wstar2, zz[2], zzq, zzt;
/* CENTRY */
int pmix_( int dyn_in_prm, float dyn_in_val, float dyn_in_val2,
float CONVECT, float CONV_CRIT, float pressure, int air_moist_prm,
float air_moist_val_dum, int sfc_moist_prm, float sfc_moist_val_dum,
float salinity, int ss_prm, float ss_val_dum, float t_air,
int sst_prm, float t_skin, float ref_ht_wind, float ref_ht_tq,
int astab, float Qnet, int warn, int flux_model, int z0_mom_prm,
int z0_TQ_prm, int stable_prm,
float *shf, float *lhf, float *tau, float *u_star,
float *t_star, float *q_star, float *z_over_L, float *wave_age,
float *dom_phs_spd, float *h_sig, float *ww_stab, float *zo_m )
/* 'main' subroutine for combined sea state and atmospheric flux model
Includes:
BVW sea state parameterization,
BVW atmospheric boundary layer (roughness length) model (optional),
Toba's 2/3 wave relation,
Monin-Obhukov atmoshpheric stability parameterization,
Godfrey and Beljaar's convective parameterization.
Options for dynamic input parameters
Most users will want to input wind speeds. However, there is a
growing need to convert scatterometer observations of stress or
equivalent neutral wind speed to winds and/or stresses.
dyn_in_prm Parameter treated as known
0 Wind speed,
1 Friction velocity,
2 Surface stress,
3 Equivalent neutral wind speed. Note 'equivalent neutral'
is NOT equivelent to 'neutral'.
Options for seastate parameterizations (or parameters):
There are three possible seastate assumptions: any one of the
following can be treated as known: wind-wave stability parameter,
phase speed, or wave age.
ss_prm Parameter treated as known
0 Wind-wave stability parameter (set to 1.0 for local
equilibrium),
1 Phase speed of the dominant waves [m/s],
2 Wave age the dominant waves (cp/ustar) [],
3 Significant Wave Height [m],
4 Significant slope [],
5 Period of the dominant waves [s],
6 Orbital velocity of dominant waves [m/s].
Options for atmospheric moisture input:
Choose the moisture parameter that is easiest for you to deal with:
air_moist_prm Parameter for moisture of air
0 Absolute humidity at the reference height of the thermometer
and humidity sensor [g vapor / g air],
1 Relative Humidity [fraction],
2 Dew point temperature [C],
3 Wet bulb temperature [C].
Options for surface moisture input:
Choose the moisture parameter that is easiest for you to deal with:
sfc_moist_prm Parameter for moisture of air
0 Absolute humidity 'at' (near) the surface [g vapor / g air],
1 Relative Humidity [fraction],
2 Dew point temperature [C],
3 Wet bulb temperature [C].
Definitions:
Input parameters (passed through call):
CON_P Convective parameter (unitless). Recammended value between
0.7 and 1.25. For details see TOGA NOTES #4.
CRIT Convergence criterion (unitless fraction).
air_moist_prm Moisture parameter index.
air_moist_val Value of the parameter corresponding to the above index.
press Atmospheric surface pressure (Pa).
s Salinity (unitless).
ss_prm Seastate parameter index.
ss_val Value of the parameter corresponding to the above index.
ta Air temperature at the reference height of the thermometer
and humidity sensor (degrees C).
u Mean wind speed (m/s) at the height (zref) of the annemometer.
u_sfc Mean surface current (m/s).
zref Height (metres) of the wind observations.
zrefq Height (metres) of the humidity observations. See zreft.
zreft Height (metres) of the temperature observations. Note: in
the current version of the code this must equal to height
of the humidity observations.
astab Option for atmospheric stability condition:
0) atmospheric stability is assumed to be neutral,
1) assumption of a neutral atmospheric stability is
not made.
theta_ang Direction of 'dynamic input' (e.g., wind or stress) relative
to the direction of propagation of the dominant waves.
Other parameters:
bstar Component of bouyancy flux (-bstar * ustar).
cmin Minimun phase speed of water waves - determined through
the phase relation [m/s].
CP Specific heat of air [J/kg].
denair density of air [kg/m^3].
denwat density of water [kg/m^3].
LV Latent heat of vaporization [J/kg].
sfcten suface tension of air-water interface [N/m].
ww Wind-wave stability parameter (unitless). Equal to one for
local equilibrium; greater than one for rising seas; less
than one (and greater than zero) for falling seas.
ww_eql Value of U(Hs)/cp for local wind-wave equilibrium (0.83).
zz wind reference height divided by momentum roughness lentgh [],
zzq moisture reference height divided by moisture roughness
lentgh [],
zzt temperature reference height divided by temperature roughness
lentgh [],
Output parameters:
cp dominant phase speed of gravity waves [m/s],
hsig significant wave height [m],
inv_monin inverse Monin-Obhukov scale length [m^-1],
lhf latent heat flux [W/m^2],
qstar scaling parameter for moisture [],
shf sensible heat flux [W/m^2],[W/m^2],
tau stress vector [N/m^2]; tau[0] is parallel the direction
of wave propagation, and tau[1] is perpendicular tau[0] in
a right handed coordinate system (with the positive
vertical axis pointing upward,
tstar scaling term for potential temperature [degrees C]
ustar friction velocity [m/s],
wa wave age, cp/ustar [],
zo_m momentum roughness lentgh (two components: parallel and
perpendicular direction of wave propagation) [],
*/
{
short fixedcp, fixedwa, neutral, fixed_hsig, fixed_uen,
fixed_ustar;
int count, z0_prm;
float viscosity;
double air_moist_val, CON_P, CRIT, cp, denwat, MS, MW, nu, qmixa,
qmixw, press, s, sfcten, sfc_moist_val,
ss_val, ta, tskin, va, ww, ww_eql, zref, zrefq,
zreft, betac_set, betag_set, betas_set, dud;
/* definitions for flux model types
0 BVW
1 Smith 1988
*/
flux_model_dum = flux_model;
smith88 = FALSE;
no_capw = FALSE;
wave_stab_intr = TRUE;
use_dh = FALSE;
if ( flux_model < 0 ) {
z0_prm = z0_mom_prm;
TQzo_prm = z0_TQ_prm;
stab_prm = stable_prm;
} else {
TQzo_prm = 0;
stab_prm = 0;
switch( flux_model ) {
case 0: /* BVW */
z0_prm = 0;
break;
case 1: /* Smith 1988 */
z0_prm = 0;
ss_prm = 2;
ss_val_dum = 43.64;
dyn_in_val2 = 0.0;
smith88 = TRUE;
no_capw = TRUE;
wave_stab_intr = FALSE;
stab_prm = 0;
break;
case 2: /* BVW without interaction between waves and atmospheric stab */
z0_prm = 0;
wave_stab_intr = FALSE;
break;
case 3: /* BVW with Smith (1988) stability parameterization */
z0_prm = 0;
stab_prm = 1;
break;
case 4: /* BVW without capillary waves */
z0_prm = 0;
no_capw = TRUE;
break;
case 5: /* BVW without sfcten in phase relations */
z0_prm = 0;
no_sfcten = TRUE;
break;
case 6: /* BVW without capillary waves and without sfcten in
the phase relation */
z0_prm = 0;
no_sfcten = TRUE;
no_capw = TRUE;
break;
case 7: /* Taylor and Yelland (2001) roughness length */
z0_prm = 2;
no_sfcten = TRUE;
no_capw = TRUE;
break;
case 8: /* Taylor and Yelland (2001) roughness length + capillary */
z0_prm = 3;
break;
case 9: /* Bourassa 2006 friction velocity with CFC zot and zoq */
z0_prm = 1;
TQzo_prm = 1;
break;
case 10: /* case 9 with displacement height estimated within the code */
z0_prm = 1;
TQzo_prm = 1;
use_dh = TRUE;
break;
case 11: /* Bourassa 2006 friction velocity with Zilitinkevich et al. zot and zoq */
z0_prm = 1;
TQzo_prm = 2;
break;
case 12: /* Bourassa 2006 friction velocity with LKB zot and zoq */
z0_prm = 1;
TQzo_prm = 3;
break;
case 13: /* Bourassa 2006 friction velocity with COARE3.0 zot and zoq */
z0_prm = 1;
TQzo_prm = 4;
break;
case 14: /* Bourassa 2006 friction velocity with wall theory for moisture */
z0_prm = 1;
TQzo_prm = 0;
break;
case 15: /* Bourassa 2006 friction velocity with CFC zot and zoq
and oil spill z0m and zero LHF */
z0_prm = 4;
TQzo_prm = 1;
break;
default:
printf( "Invalid choice of flux_model parameter: %i\n", flux_model );
exit(1);
break;
}
}
/* printf( " flux_model = %i %i %i %i\n", flux_model, z0_prm, TQzo_prm, stab_prm ); */
MS = 58.443000; /* molecular mass of salt [kg/kmol] */
MW = 18.016000; /* molecular mass of water [kg/kmol] */
betag_prime = 1.0;
betac_prime = 1.0;
if ( no_capw ) betac_prime = 0.0;
/* printf( "F1 %f\n", betac_prime ); */
betas_prime = 0.0 + (double)smith88;
betas[0] = 1.0;
betas[1] = 1.0;
ww_eql = 0.345 / KV;
CRIT = (double)CONV_CRIT;
CON_P = (double)CONVECT;
press = (double)pressure;
air_moist_val = (double)air_moist_val_dum;
sfc_moist_val = (double)sfc_moist_val_dum;
s = (double)salinity;
ss_val = (double)ss_val_dum;
ta = (double)t_air;
tskin = (double)t_skin;
zref = (double)ref_ht_wind;
/* the height of the temperature sensor must be equal to the height of
the humidity sensor. */
zrefq = (double)ref_ht_tq;
zreft = zrefq;
/* printf( "%f %lf %lf %lf\n", wind_spd, wind_spd2, CRIT, CON_P );
printf( "%lf %lf %lf\n", press, air_moist_val, sfc_moist_val );
printf( "%lf %lf %lf %lf\n", s, ss_val, ta, tskin );
printf( "%lf %lf %lf \n", zref, zreft, zrefq );
*/
/* setup the knowns and initial gueses for the seastate. */
wa[0] = 20.0; /* initial guess */
wa[1] = 28.0;
cp = 4.0; /* initial guess */
fixedcp = FALSE;
fixedwa = FALSE;
hsig[0] = 0.01; /* initial guess assumes waves exist */
hsig[1] = 0.01; /* initial guess assumes waves exist */
fixedcp = FALSE;
fixedwa = FALSE;
fixed_hsig = FALSE;
ww = 1.0;
switch( ss_prm ) {
case 0: /* the wind-wave stability parameter is specified */
ww = ss_val;
break;
case 1: /* phase speed is known */
fixedcp = TRUE;
cp = ss_val;
break;
case 2: /* wave age (ustar/cp) is known */
fixedwa = TRUE;
wa[0] = ss_val;
break;
case 3: /* significant wave height is known */
fixed_hsig = TRUE;
hsig[0] = ss_val;
break;
case 4: /* significant slope is known */
case 5: /* period of the dominant waves is known */
break;
case 6: /* Orbital velocity is known */
u_orb[0] = ss_val;
break;
default:
printf( "Invalid choice of ss_prm: %i\n", ss_prm );
exit(1);
break;
}
/* setup the atmospheric stability parameters */
monin_inv = 0.0;
switch( astab ) {
case 0:
neutral = TRUE;
break;
case 1:
neutral = FALSE;
break;
default:
printf( "Invalid atmospheric stability option: astab = %i\n", astab );
printf( "Allowable options are: 0 (neutral), and 1 (non-neutral).\n" );
exit(1);
}
betag_set = betag_prime;
betac_set = betac_prime;
betas_set = betas_prime;
/* determine the viscosity of air */
nu = 1.3260000e-5 * (1.00000 + 0.006542 * ta + 8.301000E-6 * ta * ta +
4.840000E-9 * ta * ta * ta );
/* determine the surface tension */
sfcten = 7.6100000E-2 - 1.55000E-4 * tskin + 2.77000E-5 * s * MS / MW;
/* determine the latent heat of vaporization */
LV = 4186.8 * ( 597.31 - 0.56525 * tskin );
/* use the atmospheric moisture parameter to determine the
specific humidity [kg/kg]. */
qmixa = find_q_( air_moist_prm, air_moist_val, press, ta );
if ( qmixa == -1.0 ) printf( "Invalid choice of air_moist_prm: %i\n", air_moist_prm );
/* use the surface moisture parameter to determine the
specific humidity [kg/kg]. */
qmixw = find_q_( sfc_moist_prm, sfc_moist_val, press, tskin );
if ( qmixw == -1.0 ) printf( "Invalid choice of sfc_moist_prm: %i\n", sfc_moist_prm );
if ( z0_prm == 4 ) qmixw = qmixa;
/* determine the heat capacity at constant pressure */
CP = 1004.0 * ( 1.0 + 0.9 * qmixw );
/* determine the density of the moist air */
denair = press / (R * (ta + 273.14) * (1.00000000 + 0.6100000 * qmixa ) );
/* determine the density of pure water */
denwat = (999.8396 + 18.224944 * tskin - 0.007922210 * tskin * tskin ) /
( 1.0000000 + 0.018159725 * tskin );
/* determine the density of saline water */
va = ( 12.97000 + 0.23400000 * tskin - 4.2100000E-3 * tskin * tskin +
2.8570000E-5 * tskin * tskin * tskin ) * 10.00000E-3 +
sqrt( s * denwat * 10.00000E-3 ) *
2.9820000E-3 - 4.9700000E-5 * tskin + 6.0320000E-7 * tskin * tskin;
denwat = denwat * ( 1.0000000 + s ) / ( 1.000000 + va * denwat * s / MS );
/*setup known inputs for the dynamic input (typically wind speed) */
fixed_ustar = FALSE;
fixed_uen = FALSE;
switch( dyn_in_prm ) {
case 0: /* wind speed is specified */
/* split the wind into component parallel and perpendicular to the mean
direction of wave propagation. */
u[0] = (double)dyn_in_val;
u[1] = (double)dyn_in_val2;
break;
case 1: /* friction velocity is specified */
fixed_ustar = TRUE;
/* split the wind into component parallel and perpendicular to the mean
direction of wave propagation. */
ustar[0] = (double)dyn_in_val;
ustar[1] = (double)dyn_in_val2;
if ( fabs( ustar[0] ) < 0.00001 ) ustar[0] = 0.00001;
if ( fabs( ustar[1] ) < 0.00001 ) ustar[1] = 0.00001;
break;
case 2: /* surface stress (magnitude) is known */
fixed_ustar = TRUE;
/* split the wind into component parallel and perpendicular to the mean
direction of wave propagation. */
ustar[0] = sqrt( (double)dyn_in_val / denair );
ustar[1] = sqrt( (double)dyn_in_val2 / denair );
if ( fabs( ustar[0] ) < 0.00001 ) ustar[0] = 0.00001;
if ( fabs( ustar[1] ) < 0.00001 ) ustar[1] = 0.00001;
break;
case 3: /* equivalent neutral wind speed is known */
fixed_uen = TRUE;
/* split the wind into component parallel and perpendicular to the mean
direction of wave propagation. */
u[0] = (double)dyn_in_val;
u[1] = (double)dyn_in_val2;
break;
default:
printf( "Invalid choice of dyn_in_prm: %i\n", dyn_in_prm );
exit(1);
break;
}
/* One the first pass assume than waves exist */
hsig[1] = 0.01;
zz[0] = 10000.0;
zz[1] = 10000.0;
count = solve( fixedcp, fixedwa, fixed_hsig, fixed_uen, fixed_ustar,
neutral, CON_P, cp, CRIT, denwat, nu, qmixa, qmixw, sfcten,
ss_prm, ss_val, ta, tskin, Qnet, zref, zreft, zrefq, ww, ww_eql, betag_set,
betac_set, betas_set, z0_prm, warn );
if ( count <= 1 || ustar[0] == 10.0 || ustar[1] == 10.0 || ( wa[0] == 1.08
&& betag[0] > 0.5 && flux_model != 9 ) || ( wa[0] == 250.0 && betag[0] > 0.5 &&
flux_model != 9 ) || fabs(monin_inv) == 10.0 ) {
count = -1;
if ( warn ) printf( "non-conv diags: %i %lf %lf %lf %lf %lf %lf %lf\n", count,
ustar[0], ustar[1], tstar, qstar, wa[0], monin_inv, betag[0] );
if ( !warn && ( ustar[0] == 10.0 || ustar[1] == 10.0 ) ) {
ustar[0] = 0.0; /* usually a good approximation */
ustar[1] = 0.0;
printf( "friction velocity set to zero. \n" );
}
}
if ( count > 1 || !warn ) {
dud = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
tau[0] = (float)(denair * dud * ustar[0]);
tau[1] = (float)(denair * dud * ustar[1]);
*shf = (float)(-denair * CP * dud * tstar);
*lhf = (float)(-denair * LV * dud * qstar);
u_star[0] = (float)ustar[0];
u_star[1] = (float)ustar[1];
*t_star = (float)tstar;
*q_star = (float)qstar;
*z_over_L = (float)(!neutral) * (float)( zref * monin_inv );
*wave_age = (float)wa[0];
*dom_phs_spd = (float)( betag[0] * cp_com );
*h_sig = (float)( betag[0] * hsig[0] );
*ww_stab = (float)ww;
zo_m[0] = zref / zz[0];
zo_m[1] = zref / zz[1];
} else {
/* printf( "Warning: non-convergence" ); */
tau[0] = -9999.0;
tau[1] = -9999.0;
*shf = -9999.0;
*lhf = -9999.0;
u_star[0] = (float)ustar[0];
u_star[1] = (float)ustar[1];
*t_star = (float)tstar;
*q_star = (float)qstar;
*z_over_L = (float)( zref * monin_inv );
*wave_age = (float)wa[0];
*dom_phs_spd = (float)cp_com;
*h_sig = (float)hsig[0];
*ww_stab = (float)ww;
zo_m[0] = -9999.0;
zo_m[1] = -9999.0;
}
return count;
/* end of subroutine pmix_ */
}
double find_q_( int moist_prm, double moist_val, double press, double temperature )
{
double satvp, spec_hum;
/* printf( "Bq%i %f %f %f\n", moist_prm, moist_val, press, temperature ); */
/* use the atmospheric moisture parameter to determine the
specific humidity [kg/kg]. */
switch ( moist_prm ) {
case 0: /* input is specific humidity */
spec_hum = moist_val;
break;
case 1: /* relative humidity (fraction) is used as input */
/* determine saturation vapour pressure over a smooth surface of
water, at the temperature */
satvp = (1.000700 + 3.460000E-8 * press) * 611.2100 *
exp( 17.50200 * temperature / ( 240.97000 + temperature ) );
spec_hum = moist_val * satvp * 0.6220 /
( press - 0.3780 * moist_val * satvp );
/* printf( "%f %f %f %f\n", satvp, spec_hum, temperature, press ); */
break;
case 2: /* dew point temperature is used as input */
/* determine the latent heat of vaporization */
satvp = (1.000700 + 3.460000E-8 * press) * 611.2100 *
exp( 17.50200 * moist_val / ( 240.97000 + moist_val ) );
spec_hum = 0.6220 * satvp / ( press - 0.3780 * satvp );
break;
case 3: /* wet bulb temperature is used as input */
spec_hum = 0.0000066 * ( 1.0 + 0.00115 * moist_val ); /* dummy use of spec_hum */
satvp = (1.000700 + 3.460000E-8 * press) * 611.2100 *
exp( 17.50200 * moist_val / ( 240.97000 + moist_val ) ) -
spec_hum * press * (temperature - moist_val );
spec_hum = 0.6220 * satvp / ( press - 0.3780 * satvp );
break;
default:
printf( "Invalid choice of moist_prm: %i\n", moist_prm );
spec_hum = -1.0;
}
return spec_hum;
/* end for subroutine find_q_ */
}
int ht_adj_( int dyn_in_prm, float dyn_in_val, float dyn_in_val2,
float CONVECT, float CONV_CRIT, float pressure, int air_moist_prm,
float air_moist_val, int sfc_moist_prm, float sfc_moist_val,
float salinity, int ss_prm, float ss_val, float t_air,
int sst_prm, float t_skin, float ref_ht_wind, float ref_ht_tq,
float z_wanted, int astab, int eqv_neut_prm, float Qnet, int warn,
int flux_model, int z0_mom_prm, int z0_TQ_prm, int stable_prm,
float *shf, float *lhf, float *tau, float *u_star,
float *t_star, float *q_star, float *z_over_L, float *wave_age,
float *dom_phs_spd, float *h_sig, float *ww_stab, float *zo_m,
float *u_at_z, float *t_at_z, float *q_at_z )
{
int astab_temp, bvw_flag, eqv_neut;
float alpha, CP, Dalton, denair, nu, phi_u, phi_Q, phi_T, q_sfc, renewalTime,
Stanton, ustar_mag, viscosity, zoN[2], zo_mag;
/*
printf( "ht_adj.c: %f %f %f %f %f %i %f %i %f %f %i %f %f %f %f %f %i %i %i %f\n",
dyn_in_val, dyn_in_val2,
CONVECT, CONV_CRIT, pressure, air_moist_prm, air_moist_val,
sfc_moist_prm, sfc_moist_val, salinity, ss_prm, ss_val, t_air,
t_skin, ref_ht_wind, ref_ht_tq, z_wanted, astab, warn, eqv_neut);
*/
TQzo_prm = 0;
if ( flux_model == 9 || flux_model == 10 ) TQzo_prm = 1;
if ( eqv_neut_prm == 2 ) {
eqv_neut = FALSE;
} else if ( eqv_neut_prm == 1 ) {
eqv_neut = TRUE;
} else {
eqv_neut = FALSE;
/* calculate neutral roughness length */
astab_temp = 0;
bvw_flag = pmix_( dyn_in_prm, dyn_in_val, dyn_in_val2, CONVECT, CONV_CRIT,
pressure, air_moist_prm, air_moist_val, sfc_moist_prm, sfc_moist_val,
salinity, ss_prm, ss_val, t_air, sst_prm, t_skin, ref_ht_wind, ref_ht_tq,
astab, Qnet, warn, flux_model, z0_mom_prm, z0_TQ_prm, stable_prm,
shf, lhf, tau, u_star, t_star, q_star, z_over_L, wave_age,
dom_phs_spd, h_sig, ww_stab, zo_m );
/* record neutral roughness length */
zoN[0] = zo_m[0];
zoN[1] = zo_m[1];
}
bvw_flag = pmix_( dyn_in_prm, dyn_in_val, dyn_in_val2, CONVECT, CONV_CRIT,
pressure, air_moist_prm, air_moist_val, sfc_moist_prm, sfc_moist_val,
salinity, ss_prm, ss_val, t_air, sst_prm, t_skin, ref_ht_wind, ref_ht_tq,
astab, Qnet, warn, flux_model, z0_mom_prm, z0_TQ_prm, stable_prm,
shf, lhf, tau, u_star, t_star, q_star, z_over_L, wave_age,
dom_phs_spd, h_sig, ww_stab, zo_m );
/* LV = 4186.8 * ( 597.31 - 0.56525 * t_skin ); */
q_sfc = find_q_( sfc_moist_prm, sfc_moist_val, pressure, t_skin );
/* determine the heat capacity at constant pressure */
CP = 1004.0 * ( 1.0 + 0.9 * q_sfc );
alpha = 1.0/(t_skin+273.15);
/* correct to the height of z_wanted */
if ( eqv_neut_prm == 0 ) {
phi_u = phi_u_f( z_wanted, z_wanted / zo_m[0] );
nu = 1.3260000e-5 * (1.00000 + 0.006542 * t_skin +
8.301000E-6 * t_skin * t_skin + 4.840000E-9 * t_skin * t_skin * t_skin );
ustar_mag = sqrt( ustar[0]*ustar[0] + ustar[1] * ustar[1] );
zo_mag = sqrt( zo_m[0]*zo_m[0] + zo_m[1] * zo_m[1] );
zo_mag = zo_m[0]; /* test option */
viscosity = 1.4e-5;
denair = pressure / (R * ( t_skin + 273.14) * (1.0 + 0.61 * q_sfc ) );
renewalTime = getRenewalTime( zo_mag, ustar_mag, denair, CP, alpha,
Qnet, viscosity );
Stanton = getStanton( ustar_mag, renewalTime);
Dalton = getDalton( ustar_mag, renewalTime);
/* zzt = z_wanted * z_o_zt( ref_ht_wind, ref_ht_tq, nu, zo_mag, Stanton )
/ ref_ht_tq;
zzq = z_wanted * z_o_zq( ref_ht_wind, ref_ht_tq, nu, zo_mag, Dalton ) /
ref_ht_tq; */
zzt = zzt * z_wanted / ref_ht_tq;
zzq = zzq * z_wanted / ref_ht_tq;
phi_T = phi_T_f( z_wanted, zzt );
phi_Q = phi_Q_f( z_wanted, zzq );
} else if ( eqv_neut_prm == 1 ) { /* Tang and Liu eqv neut winds */
phi_u = 0.0;
phi_T = 0.0;
phi_Q = 0.0;
} else if ( eqv_neut_prm == 2 ) {
/* if option was selected, calculate equivalent friction velocities as
defined by Geernaert and Katsaros (1986) */
ustar[0] = pow( dyn_in_val / ustar[0] + phi_u / KV -
log( zoN[0] / zo_m[0] ) / KV, -1.0 ) * dyn_in_val;
ustar[1] = pow( dyn_in_val2 / ustar[1] + phi_u / KV -
log( zoN[1] / zo_m[1] ) / KV, -1.0 ) * dyn_in_val2;
phi_u = 0.0;
}
*z_over_L = *z_over_L * z_wanted / ref_ht_wind;
*u_at_z = sqrt( pow(ustar[0] * ( log( z_wanted / zo_m[0] ) + phi_u ) / KV, 2.0) +
pow(ustar[1] * ( log( z_wanted / zo_m[1] ) + phi_u ) / KV, 2.0) );
*t_at_z = t_skin + Prt * tstar * ( log( zzt ) + phi_T ) / KV;
*q_at_z = q_sfc + Sc * qstar * ( log( zzq ) + phi_Q ) / KV;
/*
printf( "ht_adj.c: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
*shf, *lhf, tau[0], tau[1], u_star[0], u_star[1], *t_star, *q_star,
*z_over_L, *wave_age, *dom_phs_spd, *h_sig, *ww_stab, *zo_m, *u_at_z,
*t_at_z, *q_at_z, phi_u, phi_T, phi_Q );
printf( "%f %f %f %f %f %f %f\n",
z_wanted, zzt, zzq, phi_T, phi_Q, tstar, qstar );
*/
/* printf( "%i %f %f %f\n", flux_model, renewalTime, Stanton, Dalton);
printf( "%f, %f, %f %f\n", 5.0,-1.0/(Prt*Stanton),5.0,-1.0/(Sc*Dalton));
*/
/*
printf("In solve: zzt = %f; zzq = %f; phi_Q = %f; phi_T = %f\n", zzt, zzq, phi_Qx, phi_Tx );
printf( "In ht_adj: zzt = %f; zzq = %f; phi_Q = %f; phi_T = %f\n", zzt*ref_ht_tq/z_wanted, zzq*ref_ht_tq/z_wanted, phi_Q, phi_T );
printf( "test2: %f %f %f\n", betas[0], betac[0], betag[0] );
*/
return( bvw_flag );
}
/* ENDCENTRY */
double f_bstar( short neutral, double qmixa, double ta, double minimum_value )
{
double bstar;
bstar = (double)(!neutral) * G * ( ( 1.00 + 0.610 * qmixa ) * tstar + 0.6100 *
( ta + 273.15 ) * qstar ) / ( ta + 273.15 * ( 1.0 + 0.610 * qmixa ) );
if ( fabs( bstar ) < minimum_value ) {
if ( bstar == 0.0 )
bstar = minimum_value;
else
bstar = minimum_value * bstar / fabs( bstar );
}
return bstar;
}
int solve( short fixedcp, short fixedwa, short fixed_hsig,
short fixed_uen, short fixed_ustar, short neutral,
double CON_P, double cp, double CRIT, double denwat, double nu,
double qmixa, double qmixw, double sfcten, int ss_prm, double ss_val,
double ta, double tskin, float Qnet, double zref, double zreft, double zrefq,
double ww, double ww_eql, double betag_set, double betac_set,
double betas_set, int z0_prm, int warn )
{
short done, first, hsig_known, unreasonable, trouble;
int count, count_Lu, count_tqL, i;
float alpha, Dalton, renewalTime, Stanton, viscosity, zo_mag;
double bstar, delta_q, delta_theta,
max, min, monin_inv_old, monin_inv_old2, phi_T, phi_Q,
qstar_old, ta_ept, tstar_old, tskin_ept,
ustar_mag, ustar_mag_old, ustar_old[2], zeta, wa_old[2], zz_old[2];
/* convert the temperatures to equivalent potential temperature */
/* estimated through the dry static energy */
tskin_ept = tskin;
ta_ept = ta + G * zreft / CP;
betag_prime = betag_set;
betag[0] = betag_prime;
betag[1] = betag_prime;
betac_prime = betac_set;
betac[0] = betac_prime;
betac[1] = betac_prime;
betas_prime = betas_set;
betas[0] = betas_prime;
betas[1] = betas_prime;
/* intial guess for value of friction velocity (ustar) */
/* component parallel direction of wave motion */
ustar[0] = (double)fixed_ustar * ustar[0] +
(double)!fixed_ustar * 0.003;
/* component perpendicular direction of wave motion */
ustar[1] = (double)fixed_ustar * ustar[1] +
(double)!fixed_ustar * 0.003;
tstar = 0.0003 * (ta - tskin);
if ( fabs( tstar ) < 0.001 ) tstar = -0.001;
qstar = 0.0003 * (qmixa - qmixw);
if ( fabs( qstar ) < 0.0001 ) qstar = -0.0001;
ustar_old[0] = ustar[0] / 2.00;
ustar_old[1] = ustar[1] / 2.00;
tstar_old = tstar / 2.00;
qstar_old = qstar / 2.00;
cmin = sqrt( 2.0 * sqrt( G * sfcten / denwat ) );
u_sfc[0] = 0.0;
u_sfc[1] = 0.0;
if ( ss_prm != 6 ) {
u_orb[0] = 0.0;
} else {
u_orb[0] = ss_val;
}
u_orb[1] = 0.0;
/*initially assume a neutral atmospheric stability (monin_inv = 0) */
monin_inv = 0.0;
unreasonable = TRUE;
while ( unreasonable == TRUE ) {
/* iterate over atmospheric stability and friction velocity */
monin_inv_old2 = -0.001;
ustar_mag_old = 0.0;
ustar_mag = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
count_Lu = 0;
while ( ( fabs( monin_inv - monin_inv_old2 ) > CRIT*fabs( monin_inv_old2 ) ||
fabs( ustar_mag - ustar_mag_old ) > CRIT * fabs( ustar_mag_old ) ||
count_Lu < 2 ) && count_Lu < 30 && !( fixed_ustar && count_Lu >= 1 ) ) {
count_Lu++;
/* determine roughness lengths for heat and moisture profiles */
viscosity = 1.4e-5;
alpha = 1.0/(ta+273.15);
zo_mag = zref / sqrt( zz[0] * zz[0] + zz[1] * zz[1] );
zo_mag = zref / zz[0]; /* test option */
renewalTime = getRenewalTime( zo_mag, ustar_mag, denair, CP, alpha,
Qnet, viscosity );
Stanton = getStanton( ustar_mag, renewalTime);
Dalton = getDalton( ustar_mag, renewalTime);
zzt = z_o_zt( zref, zreft, nu, zo_mag, Stanton );
zzq = z_o_zq( zref, zrefq, nu, zo_mag, Dalton );
ustar_mag_old = ustar_mag;
monin_inv_old2 = monin_inv;
/* iterate over solutions to tstar, qstar, and monin_inv */
tstar_old = 0.5 * tstar;
qstar_old = 0.5 * qstar;
monin_inv_old = monin_inv;
count_tqL = 0;
while ( ( fabs( tstar - tstar_old ) > CRIT * fabs( tstar_old ) ||
fabs( qstar - qstar_old ) > CRIT * fabs( qstar_old ) ||
fabs( monin_inv - monin_inv_old ) > 0.1*CRIT * fabs( monin_inv_old )
|| count_tqL < 3 ) && count_tqL < 40 ) {
count_tqL++;
tstar_old = tstar;
qstar_old = qstar;
monin_inv_old = monin_inv;
/* determine tstar (assuming the current L is correct) */
/* printf( "test2a: %f %f %f\n", (float)tstar, (float)qstar, (float)monin_inv );
printf( "test2c: %f %f %f %f %f\n", (float)delta_theta, (float)KV, (float)( log( zzt+1.0 ), (float)(!neutral) * phi_T ), (float)Prt ); */
if ( ta_ept == tskin_ept )
tstar = 0.0;
else {
/* determine the eqv. pot. temperature stability term */
phi_T = phi_T_f( zreft, zzt );
delta_theta = ta_ept - tskin_ept;
if ( fabs( delta_theta ) < 0.001 ) delta_theta = -0.001;
tstar = KV * delta_theta / ( log( zzt+1.0 ) +
(double)(!neutral) * phi_T ) / Prt;
}
/* determine qstar (assuming current L is correct */
/* printf( "test2b: %f %f %f\n", (float)tstar, (float)qstar, (float)monin_inv ); */
if ( qmixa == qmixw )
qstar = 0;
else {
/* determine the moisture stability term */
phi_Q = phi_Q_f( zreft, zzq );
delta_q = qmixa - qmixw;
if ( fabs( delta_q ) < 0.0001 ) delta_q = -0.0001;
qstar = KV * ( delta_q ) / ( log( zzq+1.0 ) +
(double)(!neutral) * phi_Q ) / Sc;
/* printf( "%f %f %f\n", delta_q, zzq, phi_Q ); */
}
/* determine the bouyancy flux */
bstar = f_bstar( neutral, qmixa, ta, 0.000001 );
/*!!!!This assumes temperature and humidity are measured at THE SAME HEIGHT */
/* determine the inverse Monin-Obhukov scale length: 1/L */
monin_inv = KV * bstar / ( ustar[0]*ustar[0] + ustar[1] * ustar[1] );
/* insure that 1/L is within reasonable bounds */
max = 1.0;
min = 0.0001;
if ( fabs(monin_inv) < min )
monin_inv = min * monin_inv / fabs( monin_inv );
if ( fabs( monin_inv ) > max )
monin_inv = max * monin_inv / fabs( monin_inv );
} /* end of iteration on tstar, qstar, and monin_inv */
/* printf( "test3: %f %f %f\n", (float)tstar, (float)qstar, (float)monin_inv ); */
done = FALSE;
count = 0;
while ( ( done == FALSE || count < 2 ) && count < 100 ) {
count++;
done = TRUE;
wa_old[0] = wa[0];
wa_old[1] = wa[1];
zz_old[0] = zz[0];
zz_old[1] = zz[1];
first = TRUE;
/* iteratively solve for u* for non-neutral conditions */
if ( !fixed_ustar ) {
find_ustar( first, fixedcp, fixedwa, fixed_hsig, fixed_uen,
neutral, count, CON_P, cp, CRIT, denwat, nu, qmixa, sfcten, ta,
ww, ww_eql, ss_prm, ss_val, zref, zreft, z0_prm );
/* printf( "test A %f %f %f %f\n", (float)ustar[0], (float)ustar[1], (float)tstar, (float)qstar ); */
} else {
i = 0;
if ( fixed_hsig && i == 0 ) hsig_known = TRUE;
else hsig_known = FALSE;
trouble = z_o_z0( fixedcp, fixedwa, neutral, hsig_known,
denwat, nu, sfcten, zref, cp, ww, ww_eql,
ss_prm, ss_val, count, i, z0_prm );
i = 1;
hsig_known = FALSE;
trouble = z_o_z0( fixedcp, fixedwa, neutral, hsig_known,
denwat, nu, sfcten, zref, cp, ww, ww_eql,
ss_prm, ss_val, count, i, z0_prm );
}
if ( ( fabs( ( wa[0] - wa_old[0]) / wa[0] ) > CRIT ||
fabs( ( wa[1] - wa_old[1] ) / wa[1] ) > CRIT ||
fabs( ( zz[0] - zz_old[0] ) / zz[0] ) > CRIT ||
fabs( ( zz[1] - zz_old[1] ) / zz[1] ) > CRIT ) && count < 60 )
done = FALSE;
}
ustar_mag = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
} /* end of iteration over friction velocity and atmospheric stability */
/* determine if the phase speed of the dominant waves is physically
reasonable. If the phase speed is less than the physical minimum,
then assume that no capillary waves are pressent. After the second
pass, if phase speed is still unacceptable then assume the surface
is smooth. */
/* if ( smith88 && monin_inv > 0 ) printf( "zo = %f %f\n", zref / zz[0], u[0] );
if ( smith88 && monin_inv > 0 ) printf( "1 %f %f %f %f %f %f %f\n", u[0], betas_prime, betac_prime, betag_prime, betas[0], betac[0], betag[0]); */
if ( fixedwa || fixedcp || ( betag[0] != 0.0 &&
wa[0] * fabs( ustar[0] ) > cmin ) || ( betag[1] != 0.0 &&
wa[1] * fabs( ustar[1] ) > cmin ) || betas_prime == 1.0 ) {
unreasonable = FALSE;
}
else
{
if ( betag[0] != 0.0 && betag[1] != 0.0 && betac_prime != 0.0 ) {
if ( fabs( u[0] - u_sfc[0] ) > fabs( u[1] - u_sfc[1] ) ) {
betag[1] = 0.0;
} else {
betag[0] = 0.0;
}
} else if ( betac_prime > 0.0 ) {
betag[0] = 1.0;
betag[1] = 1.0;
betac_prime = 0.0;
} else {
if ( betag[0] > 0.0 && wa[0] * fabs( ustar[0] ) < cmin ) {
betag[0] = 0.0 + (double)smith88;
}
if ( betag[1] > 0.0 && wa[1] * fabs( ustar[1] ) < cmin ) {
betag[1] = 0.0 + (double)smith88;
}
if ( betag[0] == 0.0 && betag[1] == 0.0 ) {
betag_prime = 0.0 + (double)smith88;
betas_prime = 1.0;
}
}
}
/* if ( smith88 && monin_inv > 0 ) printf( "2 %f %f %f %f %f %f %f\n", u[0], betas_prime, betac_prime, betag_prime, betas[0], betac[0], betag[0]); */
}
return count_Lu + fixed_ustar;
/* end of subroutine solve */
}
double f_ustar( short fixed_uen, short neutral, double bstar, double CON_P,
double phi_u, double delta_u, double ww, double zref,
double max_value, int i )
{
/* The effect of the surface current is considered in "delta_u":
delta_u = U(zref) - Usfc
This slightly reduces the speed of the algorithm. */
double d_over_zo, delta_u2, hzzg, ustar_test, wstar2;
wstar2 = pow( fabs( sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] )
* bstar ) * H, 0.6666667);
if ( monin_inv >= 0.0 ) wstar2 = 0.0;
hzzg = zz[i]; /* in this case hzzg is not related to wave height - it is a dummy variable */
if ( hzzg < 0.00001 ) hzzg = 0.00001;
if ( fabs( delta_u ) < 0.000001 ) {
ustar_test = 0.0;
} else if ( log( hzzg + 1.0 ) + (double)!fixed_uen * (double)!neutral
* phi_u > 0.01 ) {
delta_u2 = delta_u * delta_u + CON_P * wstar2;
if ( fabs( delta_u2 ) < 0.000001 ) delta_u2 = 0.000001;
delta_u = sqrt( delta_u2 ) * delta_u / fabs( delta_u );
/* calculate displacement height, multiplied by zero or one depending on
whether or not it is used. */
d_over_zo = (float)use_dh * hzzg / zref * hsig[0] * 0.8;
ustar_test = KV * delta_u / ( log( hzzg - d_over_zo + 1.0 ) +
(double)!fixed_uen * (double)!neutral * phi_u );
/* eliminate the problem of excessive stress: limit friction velocity */
if ( pow( KV / log( ( 10.0 - d_over_zo ) * zz[i] / zref + 1.0), 2.0 ) > 0.01 )
ustar_test = 0.01 * pow( log( ( 10.0 - d_over_zo ) * zz[i] / zref + 1.0) / KV, 2.0 );
if ( ustar_test > 10.0 ) ustar_test = 10.0;
}
else {
ustar_test = 10.0;
}
/* ustar is not returned because it is a global variable */
return ustar_test;
}
double phi_u_f( double zref, double zzd )
{
double a_bh, b_bh, c_bh, d_bh, phi_u, zeta, zeta0;
a_bh = 0.7;
b_bh = 0.75;
c_bh = 5.0;
d_bh = 0.35;
/* unstable case */
if ( monin_inv < 0 ) {
switch( stab_prm ) {
case 0: /* BVW */
/* find zeta( z / L ) and zeta( zo / L ) */
zeta = pow( 1.000 - 15.000 * zref * monin_inv, 0.2500 );
zeta0 = pow( 1.000 - 15.000 * zref / zzd * monin_inv, 0.2500 );
phi_u = log( (zeta0 * zeta0 + 1.0 ) * pow( zeta0 + 1.0, 2.0 ) /
( ( zeta * zeta + 1.0 ) * pow( zeta + 1.0, 2.0 ) ) ) +
2.000000 * ( atan(zeta) - atan(zeta0) );
break;
case 1: /* Smith 88 parameterization */
a_bh = sqrt( sqrt( 1.0 - 16.0 * zref * monin_inv ) );
phi_u = -2.0 * log( 0.5 * (1.0 + a_bh) ) -
log( 0.5 * (1.0 + a_bh*a_bh) ) + 2.0 *
atan( a_bh ) - 3.14159 / 2.0;
break;
}
}
else {
/* stable case */
switch( stab_prm ) {
case 0: /* BVW */
phi_u = a_bh * zref * monin_inv +
b_bh * ( zref * monin_inv - c_bh / d_bh ) *
exp(-d_bh * zref * monin_inv ) + b_bh * c_bh / d_bh;
break;
case 1: /* Smith 88 parameterization */
phi_u = 5.0 * zref * monin_inv;
break;
}
}
return phi_u;
/* end of subroutine phi_u_f */
}
double phi_Q_f( double zrefq, double zzq )
{
double a_bh, b_bh, c_bh, d_bh, phi_Q, qlamda, q0lamda;
a_bh = 1.0;
b_bh = 0.667;
c_bh = 5.0;
d_bh = 0.35;
if ( monin_inv < 0 ) {
switch( stab_prm ) {
case 0: /* BVW */
/* find zeta( z / L ) and zeta( zo / L ) */
qlamda = sqrt( 1.000 - 9.000 * zrefq * monin_inv );
q0lamda = sqrt( 1.000 - 9.000 * zrefq / zzq * monin_inv );
phi_Q = 2.0000 * log( ( q0lamda + 1 ) / ( qlamda + 1 ) );
break;
case 1: /* Smith 88 parameterization */
a_bh = sqrt( 1.0 - 16.0 * zrefq * monin_inv );
phi_Q = -2.0 * log( ( 1.0 + a_bh ) / 2.0 );
break;
}
}
else {
/* stable case */
switch( stab_prm ) {
case 0: /* BVW */
phi_Q = pow( 1.0 + 0.6667 * a_bh * zrefq * monin_inv, 1.5 ) +
b_bh * ( zrefq * monin_inv - c_bh / d_bh ) *
exp(-d_bh * zrefq * monin_inv ) + b_bh * c_bh / d_bh - 1.0;
break;
case 1: /* Smith 88 parameterization */
phi_Q = 5.0 * zrefq * monin_inv;
break;
}
}
return phi_Q;
/* end of subroutine phi_Q_f */
}
double phi_T_f( double zreft, double zzt )
{
double a_bh, b_bh, c_bh, d_bh, phi_T, tlamda, t0lamda;
a_bh = 1.0;
b_bh = 0.667;
c_bh = 5.0;
d_bh = 0.35;
if ( monin_inv < 0 ) {
switch( stab_prm ) {
case 0: /* BVW */
/* find zeta( z / L ) and zeta( zo / L ) */
tlamda = pow( 1.000 - 9.000 * zreft * monin_inv, 0.500);
t0lamda = pow( 1.000 - 9.000 * zreft / zzt * monin_inv, 0.500);
phi_T = 2.0000 * log( ( t0lamda + 1 ) / ( tlamda + 1 ) );
break;
case 1: /* Smith 88 parameterization */
a_bh = sqrt( 1.0 - 16.0 * zreft * monin_inv );
phi_T = -2.0 * log( ( 1.0 + a_bh ) / 2.0 );
break;
}
}
else {
/* stable case */
switch( stab_prm ) {
case 0: /* BVW */
phi_T = pow( 1.0 + 0.6667 * a_bh * zreft * monin_inv, 1.5 ) +
b_bh * ( zreft * monin_inv - c_bh / d_bh ) *
exp(-d_bh * zreft * monin_inv ) + b_bh * c_bh / d_bh - 1.0;
break;
case 1: /* Smith 88 parameterization */
phi_T = 5.0 * zreft * monin_inv;
break;
}
}
return phi_T;
/* end of subroutine phi_T_f */
}
double z_o_zq( double zref, double zrefq, double nu, float zo_mag,
float Dalton )
{
double betas_temp, const1, const2, Rr, ustar_mag;
/* zzq is zref / z0Q */
ustar_mag = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
if ( TQzo_prm == 0 ) { /* wall theory */
if ( betas[0] == 0.0 ) betas_temp = 1.0;
else betas_temp = betas[0];
zzq = zrefq / ( betas_temp * 0.620 * nu / fabs( ustar_mag ) );
} else if ( TQzo_prm == 1 ) { /* CFC model */
zzq = zrefq / ( zo_mag * exp( KV * ( 5.0 - 1.0 / ( Sc * Dalton ) ) ) );
} else if ( TQzo_prm == 2 ) { /* Zilitinkevich et al. 2001 */
zzq = zrefq / ( zo_mag * exp( -4.0 * sqrt( zo_mag * ustar_mag / nu ) ) );
} else if ( TQzo_prm == 3 ) { /* LKB 1979 */
Rr = zo_mag * ustar_mag / nu;
if ( Rr < 0.11 ) {
const1 = 0.177;
const2 = 0;
} else if ( Rr < 0.825 ) {
const1 = 1.376;
const2 = 0.929;
} else if ( Rr < 3. ) {
const1 = 1.026;
const2 = -0.599;
} else if ( Rr < 10.0 ) {
const1 = 1.625;
const2 = -1.018;
} else if ( Rr < 30.0 ) {
const1 = 4.661;
const2 = -1.475;
} else if ( Rr < 100.0 ) {
const1 = 34.904;
const2 = -2.067;
} else {
printf( "LKB zoq calculation failed due to large Rr: %f\n", Rr);
/* printf( "LKB z0 = %f %f %f; LKB ustar = %f; LKB nu = %f\n", zo_mag, zz[0], zz[1], ustar_mag, nu); */
const1 = 34.904;
const2 = -2.067;
}
zzq = zrefq * ustar_mag * const1 * pow( Rr, const2 ) / nu;
} else if ( TQzo_prm == 4 ) { /* COARE3.0 */
zzq = 0.000055 * pow( zo_mag * ustar_mag / nu, -0.6 ); /* zo_q */
if ( zzq < 0.00011 ) zzq = 0.00011;
zzq = zrefq / zzq;
} else if ( TQzo_prm == 5 ) { /* Modified CFC model */
zzq = zrefq / ( zo_mag * exp( KV * ( 0.5 - 1.0 / ( Sc * Dalton ) ) ) );
} else {
printf( "Invalid choice of TQzo_prm: %i\n", TQzo_prm );
exit(0);
}
return zzq;
/* end of subroutine z_o_zq */
}
double z_o_zt( double zref, double zreft, double nu, float zo_mag,
float Stanton )
{
double betas_temp, const1, const2, Rr, ustar_mag;
ustar_mag = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
if ( TQzo_prm == 0 ) {
if ( betas[0] == 0.0 ) betas_temp = 1.0;
else betas_temp = betas[0];
/* determine zreft / zT */
zzt = zreft / ( betas_temp * 0.40 * nu / fabs( ustar_mag ) );
} else if ( TQzo_prm == 1 ) {
zzt = zreft / ( zo_mag * exp( KV * ( 5.0 - 1.0 / ( Prt*Stanton ) ) ) );
} else if ( TQzo_prm == 2 ) { /* Zilitinkevich et al. 2001 */
zzt = zreft / ( zo_mag * exp( -4.0 * sqrt( zo_mag * ustar_mag / nu ) ) );
} else if ( TQzo_prm == 3 ) { /* LKB 1979 */
Rr = zo_mag * ustar_mag / nu;
if ( Rr < 0.11 ) {
const1 = 0.292;
const2 = 0;
} else if ( Rr < 0.825 ) {
const1 = 1.808;
const2 = 0.826;
} else if ( Rr < 3. ) {
const1 = 1.393;
const2 = -0.528;
} else if ( Rr < 10.0 ) {
const1 = 1.956;
const2 = -0.870;
} else if ( Rr < 30.0 ) {
const1 = 4.994;
const2 = -1.297;
} else if ( Rr < 100.0 ) {
const1 = 30.790;
const2 = -1.845;
} else {
printf( "LKB zoq calculation failed due to large Rr: %f\n", Rr);
const1 = 30.790;
const2 = -1.845;
}
zzt = zreft * ustar_mag * const1 * pow( Rr, const2 ) / nu;
} else if ( TQzo_prm == 4 ) { /* COARE3.0 */
zzt = 0.000055 * pow( zo_mag * ustar_mag / nu, -0.6 ); /* zo t */
if ( zzt < 0.00011 ) zzt = 0.00011;
zzt = zreft / zzt;
} else {
printf( "Invalid choice of TQzo_prm: %i\n", TQzo_prm );
exit(0);
}
return zzt;
/* end of subroutine z_o_zt */
}
int z_o_z0( short fixedcp, short fixedwa, short neutral, short hsig_known,
double denwat, double nu, double sfcten, double zref,
double cp, double ww, double ww_eql, int ss_prm, double ss_val, int count,
int i, int z0_prm )
{
short period_known;
int test;
double A_Charnock, B_Toba, dum, hzzg, period[2], p1, p2, p3, p4,
oil_mod, phi_u, scale, sig_slope, ueff, Uorb, ustar_temp, wave_len, z0;
A_Charnock = 0.035;
oil_mod = 0.25; /* crude guess for oil spills - not natural oil */
/* A_Charnock = 0.48/ss_val; */
test = 0;
ustar_temp = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
cmin = sqrt( 2.0 * sqrt( G * sfcten / denwat ) );
ueff = sqrt( ( u[i] - u_sfc[i] + wstar2 ) * ( u[i] - u_sfc[i] ) + wstar2 );
/* For the flux models (i.e., only BVW) that calculate sea state */
if ( z0_prm == 0 || z0_prm == 3 ) {
B_Toba = 0.0602;
/* the sea state parameter is significant slope */
if ( ss_prm == 4 && i == 0 ) {
sig_slope = ss_val;
cp = 2.0 * 3.14159 * fabs( ustar[i] ) * B_Toba * B_Toba / sig_slope / sig_slope;
wa[i] = cp / fabs( ustar[i] );
if ( wa[i] > 12000.0 ) wa[i] = 12000.0;
/* printf( "%lf %lf %lf\n", cp, wa[i], ustar[i] ); */
ueff = sqrt( u[i] * u[i] + wstar2 );
z0 = betas[i] * 0.1100 * nu / fabs( ustar[i] ) +
betac[i] * B * sfcten / ( ustar[i] * ustar[i] * denwat ) +
betag[i] * 0.016 * ustar[i] * ustar[i] / G;
phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0 );
ww = ( log( hsig[i] / z0 + 1.0 ) + phi_u ) / ( KV * ww_eql * wa[i] );
/* determine the wave characteristics */
period[i] = pow( pow( cp * sig_slope / B_Toba, 2.0 ) / ( G * ustar[i] ),
0.333333 );
hsig[i] = B_Toba * sqrt( G * fabs( ustar[i] ) * pow( period[i], 3 ) );
wave_len = cp * period[i];
}
/* the sea state parameter is significant wave height OR period */
if ( ( hsig_known || ss_prm == 5 ) && i == 0 ) {
if ( ss_prm == 5 ) period_known = TRUE;
else period_known = FALSE;
period[i] = (double)(!period_known) *
pow( hsig[i] * hsig[i] / ( G * fabs( ustar[i] ) * B_Toba*B_Toba ), 0.333333 ) +
(double)period_known * ss_val;
hsig[i] = (double)period_known * B_Toba *
sqrt( G * fabs( ustar[i] ) * pow( period[i], 3 ) )
+ (double)(!period_known) * ss_val;
p1 = G * period[i] / ( 2.0 * 3.14159 );
p2 = pow( cmin, 4.0);
p3 = 0.5 * p2 / p1 + pow( p1, 3.0 ) / 27.0;
p4 = 0.25 * p2 * p2 / ( p1 * p1 ) + p1 * p1 * p2 / 27.0;
cp = pow( p3 + sqrt( p4 ), 0.333333 ) +
pow( p3 - sqrt( p4 ), 0.333333 ) + p1 / 3.0;
wa[i] = cp / fabs( ustar[i] );
if ( wa[i] > 12000.0 ) wa[i] = 12000.0;
/* determine the wave characteristics */
wave_len = cp * period[i];
ueff = sqrt( u[i] * u[i] + wstar2 );
if ( z0_prm == 0 ) { /* BVW */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * 0.48 / wa[i] * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 1 ) { /* Bourassa (2006) */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 2 ) { /* BVW with Taylor & Yelland zog */
zz[i] = zref / ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ), 2.0 ) ) );
} else if ( z0_prm == 3 ) { /* Taylor & Yelland */
zz[i] = zref / ( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ) );
} else if ( z0_prm == 4 ) { /* Bourassa (2006) modified for oil slick */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
oil_mod * sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 5 ) { /* Smooth surface */
zz[i] = zref / ( 0.1100 * nu / fabs( ustar[i] ) );
} else {
printf( "Unreasonable value of z0_prm: %i\n", z0_prm );
}
phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0 );
ww = ( log( hsig[i] / z0 + 1.0 ) + phi_u ) / ( KV * ww_eql * wa[i] );
/* printf( "z/z0 T: i = %i; ss_val=%f; period = %f; hsig = %f; wa = %f\n", i, ss_val, period[i], hsig[i], wa[i] );
*/
} else {
ueff = sqrt( u[i] * u[i] + wstar2 );
if ( z0_prm == 0 ) { /* BVW */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * 0.48 / wa[i] * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 1 ) { /* Bourassa (2006) */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 2 ) { /* BVW with Taylor & Yelland zog */
zz[i] = zref / ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ), 2.0 ) ) );
} else if ( z0_prm == 3 ) { /* Taylor & Yelland */
zz[i] = zref / ( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ) );
} else if ( z0_prm == 4 ) { /* Bourassa (2006) modified for oil slick */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
oil_mod * sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 5 ) { /* Smooth surface */
zz[i] = zref / ( 0.1100 * nu / fabs( ustar[i] ) );
} else {
printf( "Unreasonable value of z0_prm: %i\n", z0_prm );
}
/* determine the wave characteristics */
if ( betag[i] > 0.1 ) {
dum = ustar[i] * ustar[i] * wa[i] * wa[i];
if ( dum > cmin * cmin )
wave_len = 3.14159 * ( dum + ( sqrt( dum * dum -
(double)(!no_sfcten) * pow(cmin,4) ) ) ) / G;
else
wave_len = PI2 * sqrt( sfcten / ( G * denwat ) );
period[i] = wave_len / ( fabs( wa[i] * ustar[i] ) );
/* Toba's relation between wave height, friction velocity and period */
hsig[i] = 0.062 * sqrt( G * fabs( ustar[i] ) * pow( period[i], 3 ) );
/* limit the height of breaking waves */
if ( hsig[i] > 0.137 * wave_len ) hsig[i] = 0.137 * wave_len;
} else {
hsig[i] = 0.0;
period[i] = -9999.9;
}
}
/* printf( "z/z0 Hs: i=%i; hsig=%f; ustar=%f; wa=%f; period=%f\n",
i, hsig[i], ustar[i], wa[i], period[i] );
*/
} /* end of block estimating seastate for BVW */
/* determine beta-c */
betac[i] = betac_prime;
if ( z0_prm != 3 ) {
scale = 0.8;
hzzg = scale * hsig[i] * zz[i] / zref; /* can't be -ve */
phi_u = phi_u_f( betag[i] * scale * hsig[i], hzzg );
if ( ss_prm == 6 ) {
Uorb = u_orb[i];
} else {
Uorb = 0.0 * (double)wave_stab_intr * betag[i] * scale *
0.5 * PI2 * hsig[i] / period[i];
}
phi_u = phi_u_f( 10.0, hzzg );
ueff = ustar[i] * ( log( zz[i] + 1.0 ) + phi_u ) / KV;
if ( z0_prm == 0 ) {
if ( ueff < 1.8 ) {
betac[i] = 0.0;
betag[i] = 0.0;
} else {
betac[i] = 1.0;
betag[i] = 1.0;
}
}
if ( z0_prm == 1 ) {
if ( ueff < 1.0 ) {
betac[i] = 0.0;
betag[i] = 0.0;
} else {
betac[i] = tanh( pow( 0.4 * ( ueff -1.0 ), 3 ) ) * betac_prime;
betag[i] = tanh( pow( 0.3 * ( ueff -1.0 ), 3 ) ) * betag_prime;
}
}
if ( z0_prm == 4 ) {
if ( ueff < 7.0 ) {
betac[i] = 0.0;
} else {
betac[i] = tanh( pow( 0.4 * ( ueff -7.0 ), 3 ) ) * betac_prime;
betag[i] = tanh( pow( 0.3 * ( ueff -7.0 ), 3 ) ) * betag_prime;
}
}
if ( betac[i] < 0.1 && !no_capw ) {
betag[i] = betac[i];
}
if ( no_capw ) betac[i] = 0.0;
betas[i] = 1.0 - betac[i];
/* printf( "betac: hzzg=%f; phi_u=%f; Uorb=%f\n", hzzg, phi_u, Uorb);
*/
/* determine wave age */
phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0 );
hzzg = hsig[i] * zz[i] / zref; /* can't be -ve */
if ( hzzg < 0.00001 ) hzzg = 0.00001;
if ( ss_prm == 0 || i == 1 ) {
/* local equilibrium requires ww = 1 */
wa[i] = ( log( hzzg + 1.0 ) + (double)(!neutral) * phi_u )
/ ( (double)( (1.0-i) * ww + i ) * KV * ww_eql );
cp = wa[i] * fabs( ustar[i] );
/* printf( "test1: %i %f %i %f %f %f %f %f\n", i, hzzg, neutral,
phi_u, KV, ww, ww_eql, wa[i] ); */
}
}
/* printf( "test2: %i %f %f %f %f\n", i, ueff, betas[i], betac[i], betag[i] ); */
if ( !fixedwa ) { /* i must equal 0 to reach this point */
/* wave age based on sea state observations (other than wave age) */
if ( fabs( ustar[i] ) > 0.0 ) wa[i] = cp / fabs( ustar[i] );
else wa[i] = 999.0;
} else {
cp = wa[i] * fabs( ustar[i] );
}
/* printf( "z_o_z0: i=%i; u=%f; wa=%f; ustar=%f; hsig=%f; hzzg=%f; phi_u=%f; cp=%f; bc=%f\n",
i, u[i], wa[i], ustar[i], hsig[i], hzzg, phi_u, cp, betac[i] );
*/
/* insure that wa is not horribly unreasonable */
if ( !( fixedwa || hsig_known || ss_prm == 4 ) && wa[i] > 250.0 ) {
wa[i] = 250.0;
test = -1;
}
else test = test - test%2;
if ( !fixedwa && wa[i] < 1.08 ) {
wa[i] = 1.08;
test = test - 2;
}
else test = test - ( test%4 - test%2 );
if ( !fixedcp ) cp_com = wa[0] * fabs( ustar[0] );
else cp_com = cp;
/* printf( "betac & betag: %f %f %f %f\n", betac[0], betag[1], betac[1], betag[1] );
*/
/*
printf( "betac & betag: %f %f %f %f\n", betac[i], betag[i], betac_prime, betag_prime );
*/
/* determine zref / z0 */
if ( z0_prm == 0 ) { /* BVW */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * 0.48 / wa[i] * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
/* * exp( -KV * Uorb / fabs( ustar[i] ) */
} else if ( z0_prm == 1 ) { /* Bourassa (2006) */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 2 ) { /* BVW with Taylor & Yelland zog */
zz[i] = zref / ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ), 2.0 ) ) );
} else if ( z0_prm == 3 ) { /* Taylor & Yelland */
zz[i] = zref / ( betag[i] *
hsig[i] * 1200.0 * pow( hsig[i] / wave_len, 4.5 ) );
} else if ( z0_prm == 4 ) { /* Bourassa (2006) modified for oil slick */
zz[i] = zref / ( ( ( betas[i] * 0.1100 * nu / fabs( ustar[i] ) ) +
oil_mod * sqrt( pow( betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) *
denwat ), 2.0 ) + pow( betag[i] * A_Charnock * ustar_temp *
fabs( ustar[i] ) / G, 2.0 ) ) ) );
} else if ( z0_prm == 5 ) { /* Smooth surface */
zz[i] = zref / ( 0.1100 * nu / fabs( ustar[i] ) );
} else {
printf( "Unreasonable value of z0_prm: %i\n", z0_prm );
}
/* insure that the roughness length is not too absurd */
if ( zref > zz[i] * 25 ) {
zz[i] = zref / 25.0;
test = test - 4;
}
else test = test - ( test%8 - test%4 - test%2 );
/* printf( "z_o_z0: i=%i; zz=%f; u*=%f; wa=%f; hs=%f\n", i, zz[i], ustar[i], wa[i], hsig[i] ); */
return test; /* returns an error code */
/* end of subroutine z_o_z0 */
}
void find_ustar( short first, short fixedcp, short fixedwa, short fixed_hsig,
short fixed_uen, short neutral, int count, double CON_P,
double cp, double CRIT, double denwat, double nu, double qmixa, double sfcten,
double ta, double ww, double ww_eql, int ss_prm, double ss_val, double zref,
double zreft, int z0_prm )
{
short done, done_outer, hsig_known;
int i, max, min, test, ucount;
double bstar, phi_u, ustar_old[2], ustar_temp;
ustar_temp = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] );
ucount = 0;
/* iteratively solve for u* for non-neutral conditions */
/* determine the bouyancy flux */
done_outer = FALSE;
while ( !done_outer && ucount < 30 ) {
ucount++;
for ( i=0; i<2; i++ ) {
ustar_old[i] = ustar[i];
/* ustar_temp is ustar squared */
ustar_temp = ustar[0] * ustar[0] + ustar[1] * ustar[1];
if ( fixed_hsig && i == 0 ) hsig_known = TRUE;
else hsig_known = FALSE;
bstar = f_bstar( neutral, qmixa, ta, 0.000001 );
/* !!!! This assumes temperature and humidity are measured at THE SAME HEIGHT
*/
/* determine the inverse Monin-Obhukov scale length: 1/L */
monin_inv = KV * bstar / ustar_temp;
/* insure that 1/L is within reasonable bounds */
max = 1.0;
min = 0.0001;
if ( fabs(monin_inv) < min )
monin_inv = min * monin_inv / fabs( monin_inv );
if ( fabs( monin_inv ) > max )
monin_inv = max * monin_inv / fabs( monin_inv );
if ( fabs( u[i] - u_sfc[i] ) < 0.001 ) {
ustar[i] = 0.0;
} else {
if ( !first ) done = FALSE;
first = FALSE;
test = z_o_z0( fixedcp, fixedwa, neutral, hsig_known,
denwat, nu, sfcten, zref, cp, ww, ww_eql,
ss_prm, ss_val, count, i, z0_prm );
phi_u = phi_u_f( zref, zz[i] );
/* printf( "A %f %f %f %f\n", (float)ustar[0], (float)ustar[1], (float)tstar, (float)qstar ); */
ustar[i] = f_ustar( fixed_uen, neutral, bstar, CON_P, phi_u,
u[i] - u_sfc[i], ww, zref, 10.0, i );
}
}
if ( ustar[0] <= 0.000001 ) done_outer =
fabs( ustar[0] - ustar_old[0] ) < 0.0001;
else done_outer = fabs( ( ustar[0] - ustar_old[0] )
/ ustar[0] ) < 0.01;
if ( ustar[1] <= 0.000001 ) done_outer = done_outer &&
fabs( ustar[1] - ustar_old[1] ) < 0.0001;
else done_outer = done_outer &&
fabs( ( ustar[1] - ustar_old[1] ) / ustar[1] ) < 0.01;
}
}
float getRenewalTime( float zo_mag, float ustar_mag, float density, float cp,
float alpha, float Qnet, float viscosity )
{
float renewalTime, Rfcr, Rfo, Cshear, Cconv, shearContribution,
convContribution;
Rfcr = -2.0e-04;
if(density < 50) Rfcr = -2.0e-01;
if(density < 50) Cshear = 53.32;
else Cshear = 209; /** from Wick thesis **/
if(density < 50) Cconv = 0.80;
else Cconv = 3.13; /** from Wick thesis **/
/** get surface Richardson number **/
Rfo=alpha*G*Qnet*viscosity/(density*cp*pow(ustar_mag,4.0));
if(Rfo == fabs(Rfo)) Rfo *= -1.0;
/** now get time rate **/
shearContribution = Cshear*sqrt(viscosity*zo_mag/pow(ustar_mag, 3.0));
convContribution =
Cconv*sqrt(viscosity*density*cp/(alpha*G*fabs(Qnet)));
renewalTime = (convContribution-shearContribution)*exp(-Rfcr/Rfo);
renewalTime += shearContribution;
return(renewalTime);
}
float getSolar( float solar )
{
/** this just determines the amount of solar radiation is
deposited in the skin. Aidong: this isn't needed in new
surface skin model**/
float gamma1, gamma2;
float R2;
float z;
float solarAbs;
if(waterType == 0) return(0.0);
else if(waterType == 1) {
R2 = 0.58;
gamma1 = 0.35;
gamma2 = 23;
}
else if(waterType == 1.1) {
R2 = 0.62;
gamma1 = 0.60;
gamma2 = 20;
}
else if(waterType == 1.2) {
R2 = 0.67;
gamma1 = 1.00;
gamma2 = 17;
}
else if(waterType == 2) {
R2 = 0.77;
gamma1 = 1.50;
gamma2 = 14;
}
else if(waterType == 3) {
R2 = 0.78;
gamma1 = 1.40;
gamma2 = 7.9;
}
z = 0.001;
solarAbs = solar - solar*(R2*exp(-z/gamma1) + (1-R2)*exp(-z/gamma2));
return(solarAbs);
}
float getStanton( float ustar_mag, float renewalTime)
{
float Stanton;
Stanton = sqrt(thermalDiff/(ustar_mag * ustar_mag * renewalTime));
return( Stanton );
}
float getDalton( float ustar_mag, float renewalTime )
{
float Dalton;
Dalton = sqrt(molecularDiff/(ustar_mag*ustar_mag*renewalTime));
return(Dalton);
}
float getSkin( float zo_mag, float ustar_mag, float bulkTemp, float rl, float Hs,
float Hl, float rs)
{
float Qsolar;
float Qnet;
float roWater;
float CpWater;
float uStarWater;
float roAir;
float skinTemperature;
float kappawater;
float thermalExpansion;
float viscosity;
float renewalTime;
roWater = 1021.7;
CpWater = 4002.0;
roAir = 1.25;
kappawater = 1.4e-07; /** from Gill page 71 **/
thermalExpansion = 3196e-07; /** from Gill Table A3.1 **/
viscosity = 1e-06; /** this is kinematic viscosity (Gill pg. 75) **/
/** get solar radiation absorbed over skin **/
Qsolar = getSolar(rs);
/** Get Net Heat Loss from ocean **/
Qnet = rl + Hs + Hl + Qsolar;
/** get u* for water **/
uStarWater = sqrt(roAir/roWater)*ustar_mag;
/** now get skin temperature **/
renewalTime = getRenewalTime(zo_mag,uStarWater, roWater, CpWater,
thermalExpansion, Qnet, viscosity);
skinTemperature = sqrt(renewalTime);
skinTemperature *= Qnet/(roWater*CpWater*sqrt(kappawater));
skinTemperature = skinTemperature + bulkTemp;
return(skinTemperature);
}