/* '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 (July 1998) a reasearch associate at COAPS. 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 "bvw99.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. */ /* CENTRY */ int pmix_( int dyn_in_prm, float dyn_in_val, float wind_ang, 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, float t_skin, float ref_ht_wind, float ref_ht_tq, int astab, int warn, 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 []. 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. Relative to the mean surface current. 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. rel_wnd_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.81). 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, unanswered, fixed_hsig, fixed_uen, fixed_ustar; int count, i, ID, test; double air_moist_val, *betac, *betag, *betas, CON_P, CP, CRIT, *cp, denair, denwat, *hsig, LV, *monin_inv, MS, MW, nu, phi_u, qmixa, qmixw, *qstar, press, rel_wnd_ang, s, sfcten, sfc_moist_val, ss_val, ta, tskin, *tstar, *ustar, va, *wa, *wind_vel, ww, ww_eql, zref, zrefq, zreft, *zz, betac_set, betag_set, betas_set, dud; betac = (double *) malloc( 2 * sizeof( double ) ); betag = (double *) malloc( 2 * sizeof( double ) ); betas = (double *) malloc( 2 * sizeof( double ) ); cp = (double *) malloc( sizeof( double ) ); hsig = (double *) malloc( 2 * sizeof( double ) ); monin_inv =(double *) malloc( sizeof( double ) ); qstar = (double *) malloc( sizeof( double ) ); tstar = (double *) malloc( sizeof( double ) ); ustar = (double *) malloc( 2 * sizeof( double ) ); wa = (double *) malloc( 2 * sizeof( double ) ); wind_vel = (double *) malloc( 2 * sizeof( double ) ); zz = (double *) malloc( 2 * sizeof( double ) ); MS = 58.443000; /* molecular mass of salt [kg/kmol] */ MW = 18.016000; /* molecular mass of water [kg/kmol] */ betag_set = 1.0; betac_set = 1.0; betas_set = 0.0; betas[0] = 1.0; betas[1] = 1.0; ww_eql = 0.324 / KV; /* printf( "bvw98.c: %i %f %f %f %f %f %i %f %i %f %f %i %f %f %f %f %f %i %i\n", dyn_in_prm, dyn_in_val, wind_ang, CONVECT, CONV_CRIT, pressure, air_moist_prm, air_moist_val_dum, sfc_moist_prm, sfc_moist_val_dum, salinity, ss_prm, ss_val_dum, t_air, t_skin, ref_ht_wind, ref_ht_tq, astab, warn ); */ rel_wnd_ang = (double)wind_ang; while ( wind_ang < 0.0 ) wind_ang = wind_ang + 360.0; 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, rel_wnd_ang, 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 ); printf( "%f %lf %lf %lf %lf %lf %lf\n", wind_spd, rel_wnd_ang, DTR, cos( DTR * rel_wnd_ang ), sin( DTR * rel_wnd_ang ), u[0], u[1] ); */ /* setup the knowns and initial gueses for the seastate. */ wa[0] = 20.0; /* initial guess */ wa[1] = 28.0; /* initial guess */ *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; 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; 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); } /* 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 ); /* 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. */ wind_vel[0] = (double)dyn_in_val * cos( (double)(DTR * rel_wnd_ang) ); wind_vel[1] = (double)dyn_in_val * sin( (double)(DTR * rel_wnd_ang) ); 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 * cos( (double)(DTR * rel_wnd_ang) ); ustar[1] = (double)dyn_in_val * sin( (double)(DTR * rel_wnd_ang) ); 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 ) * cos( (double)(DTR * rel_wnd_ang) ); ustar[1] = sqrt( (double)dyn_in_val / denair ) * sin( (double)(DTR * rel_wnd_ang) ); 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. */ wind_vel[0] = (double)dyn_in_val * cos( (double)(DTR * rel_wnd_ang) ); wind_vel[1] = (double)dyn_in_val * sin( (double)(DTR * rel_wnd_ang) ); 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, ustar, qstar, tstar, monin_inv, zz, CON_P, cp, wa, hsig, CRIT, denwat, nu, qmixa, qmixw, sfcten, ss_prm, ss_val, ta, tskin, wind_vel, zref, zreft, zrefq, ww, ww_eql, betag_set, betac_set, betas_set, betac, betag, betas, warn ); if ( count <= 1 || ustar[0] == 10.0 || ustar[1] == 10.0 || wa[0] == 1.08 && betag[0] > 0.5 || ( wa[0] == 250.0 && betag[0] > 0.5 ) || fabs(*monin_inv) == 1.0 ) { count = -1; if ( !warn ) printf( "non-conv diags: %i %lf %lf %lf, %lf\n", count, ustar[0], ustar[1], wa[0], betag[0] ); } 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) ); *h_sig = (float)( betag[0] * hsig[0] ); phi_u = phi_u_f( betag[0] * hsig[0], hsig[0] / zref * zz[0], monin_inv ); ww = ( log( hsig[0] / zref * zz[0] + 1.0 ) + phi_u ) / ( KV * ww_eql * wa[0] ); *ww_stab = (float)ww; zo_m[0] = (float)( zref / zz[0] ); zo_m[1] = (float)( 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); *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; /* 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.0007 + 3.46E-8 * *press) * 611.21 * exp( 17.502 * *temperature / ( 240.97 + *temperature ) ); spec_hum = *moist_val * satvp * 0.622 / ( *press - 0.378 * *moist_val * satvp ); break; case 2: /* dew point temperature is used as input */ /* determine the latent heat of vaporization */ satvp = (1.0007 + 3.46E-8 * *press) * 611.21 * exp( 17.50200 * *moist_val / ( 240.97 + *moist_val ) ); spec_hum = 0.622 * satvp / ( *press - 0.378 * 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 wind_ang, 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, float t_skin, float ref_ht_wind, float ref_ht_tq, int astab, int warn, 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, int eqv_neut, float z_wanted, float *u_at_z, float *t_at_z, float *q_at_z ) { int bvw_flag; float LV, nu, phi_u, phi_Q, phi_T, q_sfc; double monin_inv, press, sfc_mst_val, sst, ustar, zzq, zzt; /* printf( "ht_adj_: %f %f %f %f %f %i %f %i %f %f %i %f %f %f %f %f %i %i %i %f\n", dyn_in_val, wind_ang, 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, astab, warn, eqv_neut, z_wanted); */ bvw_flag = pmix_(dyn_in_prm, dyn_in_val, wind_ang, 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, astab, warn, shf, lhf, tau, u_star, t_star, q_star, z_over_L, wave_age, dom_phs_spd, h_sig, ww_stab, zo_m); /* correct to the height of z_wanted */ if ( !eqv_neut ) { monin_inv = (double)(*z_over_L / ref_ht_wind); ustar = (double)(*u_star); phi_u = phi_u_f( z_wanted, z_wanted / zo_m[0], &monin_inv ); nu = 1.3260000e-5 * (1.00000 + 0.006542 * t_air + 8.301000E-6 * t_air * t_air + 4.840000E-9 * t_air * t_air * t_air ); zzt = z_wanted * z_o_zt( ref_ht_wind, ref_ht_tq, &ustar, nu ) / ref_ht_wind; phi_T = phi_T_f( z_wanted, &zzt, &monin_inv ); zzq = z_wanted * z_o_zq( ref_ht_wind, ref_ht_tq, &ustar, nu ) / ref_ht_wind; phi_Q = phi_Q_f( z_wanted, &zzq, &monin_inv ); } else { phi_u = 0.0; phi_T = 0.0; phi_Q = 0.0; } *z_over_L = *z_over_L * z_wanted / ref_ht_wind; LV = 4186.8 * ( 597.31 - 0.56525 * t_skin ); sfc_mst_val = (double)sfc_moist_val; press = (double)pressure; sst = (double)t_skin; q_sfc = find_q_( &sfc_moist_prm, &sfc_mst_val, &press, &sst ); *u_at_z = sqrt( pow(u_star[0] * ( log( z_wanted / zo_m[0] ) + phi_u ) / KV, 2.0) + pow(u_star[1] * ( log( z_wanted / zo_m[1] ) + phi_u ) / KV, 2.0) ); *t_at_z = t_skin + *t_star * ( log( zzt ) + phi_T ) / KVT; *q_at_z = q_sfc + *q_star * ( log( zzq ) + phi_Q ) / KVQ; /* printf( "ht_adj_: %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 ); */ return( bvw_flag ); /* end of routine ht_adj_ */ } /* ENDCENTRY */ double f_bstar( short neutral, double qmixa, double ta, double *qstar, double *tstar, double minimum_value ) { double bstar; bstar = !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 *ustar, double *qstar, double *tstar, double *monin_inv, double *zz, double CON_P, double *cp, double *wa, double *hsig, double CRIT, double denwat, double nu, double qmixa, double qmixw, double sfcten, int ss_prm, double ss_val, double ta, double tskin, double *wind_vel, double zref, double zreft, double zrefq, double ww, double ww_eql, double betag_set, double betac_set, double betas_set, double *betac, double *betag, double *betas, int warn ) { short done, first, hsig_known, unreasonable, trouble; int count, count_Lu, count_tqL, i, test[2]; double betac_prime, betag_prime, betas_prime, bstar, cmin, CP, delta_u2, delta_q, delta_theta, hzzg, max, min, monin_inv_old, monin_inv_old2, phi_u, phi_T, phi_Q, qstar_old, ta_ept, tstar_old, tskin_ept, ustar_mag, ustar_mag_old, ustar_old[2], zeta, zeta0, wa_old[2], zz_old[2], *zzq, *zzt; zzq = (double *) malloc( sizeof( double ) ); zzt = (double *) malloc( sizeof( double ) ); /* convert the temperatures to equivalent potential temperature */ /* estimated through the dry static energy */ tskin_ept = tskin; CP = 1004.0 * ( 1.0 + 0.9 * qmixw ); 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 ) ); /*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 */ *zzt = z_o_zt( zref, zreft, ustar, nu ); *zzq = z_o_zq( zref, zrefq, ustar, nu ); 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 ) > 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) */ if ( ta_ept == tskin_ept ) *tstar = 0.0; else { /* determine the eqv. pot. temperature stability term */ phi_T = phi_T_f( zreft, zzt, monin_inv ); delta_theta = ta_ept - tskin_ept; if ( fabs( delta_theta ) < 0.001 ) delta_theta = -0.001; *tstar = KVT * delta_theta / ( log( *zzt+1.0 ) + !neutral * phi_T ); } /* determine qstar (assuming current L is correct */ if ( qmixa == qmixw ) *qstar = 0; else { /* determine the moisture stability term */ phi_Q = phi_Q_f( zreft, zzq, monin_inv ); delta_q = qmixa - qmixw; if ( fabs( delta_q ) < 0.0001 ) delta_q = -0.0001; *qstar = KVQ * ( delta_q ) / ( log( *zzq+1.0 ) + !neutral * phi_Q ); } /* determine the bouyancy flux */ bstar = f_bstar( neutral, qmixa, ta, qstar, tstar, 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 */ 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, ustar, qstar, tstar, monin_inv, zz, cp, wa, hsig, CRIT, denwat, nu, qmixa, sfcten, ta, wind_vel, betac_prime, betag_prime, betas_prime, betac, betag, betas, ww, ww_eql, ss_prm, ss_val, zref, zreft ); 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, betac_prime, betag_prime, betas_prime, betac, betag, betas, wind_vel, cp, wa, hsig, ww, ww_eql, ss_prm, ss_val, ustar, monin_inv, zz, count, i ); i = 1; hsig_known = FALSE; trouble = z_o_z0( fixedcp, fixedwa, neutral, hsig_known, denwat, nu, sfcten, zref, betac_prime, betag_prime, betas_prime, betac, betag, betas, wind_vel, cp, wa, hsig, ww, ww_eql, ss_prm, ss_val, ustar, monin_inv, zz, count, i ); } 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 ( fixedwa || fixedcp || ( betag[0] != 0.0 && wa[0] * fabs( ustar[0] ) > cmin ) || ( betag[1] != 0.0 && wa[1] * fabs( ustar[1] ) > cmin ) || ( betag_prime == 0.0 && betac_prime == 0.0 ) ) { unreasonable = FALSE; } else { if ( betag[0] != 0.0 && betag[1] != 0.0 && betac_prime != 0.0 ) { if ( fabs( wind_vel[0] ) > fabs( wind_vel[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; } if ( betag[1] > 0.0 && wa[1] * fabs( ustar[1] ) < cmin ) { betag[1] = 0.0; } if ( betag[0] == 0.0 && betag[1] == 0.0 ) { betag_prime = 0.0; betas_prime = 1.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 *ustar, double *monin_inv, double *zz, double max_value, int i ) { /* The effect of the surface current is considered in "delta_u": delta_u = U(zref) - Usfc This slightly increases the speed of the algorithm. */ double 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 ); ustar_test = KV * delta_u / ( log( hzzg + 1.0 ) + (double)!fixed_uen * (double)!neutral * phi_u ); /* eliminate the problem of excessive stress: limit friction velocity */ 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 *monin_inv ) { 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 ) { /* 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) ); } else /* stable case */ 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; return phi_u; /* end of subroutine phi_u_f */ } double phi_Q_f( double zrefq, double *zzq, double *monin_inv ) { 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 ) { /* 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 ) ); } else /* stable case */ 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; return phi_Q; /* end of subroutine phi_T_f */ } double phi_T_f( double zreft, double *zzt, double *monin_inv ) { 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 ) { /* 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 ) ); } else /* stable case */ 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; return phi_T; /* end of subroutine phi_T_f */ } double z_o_zq( double zref, double zrefq, double *ustar, double nu ) { double ustar_test, zzq; ustar_test = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] ); /* determine zref / z0Q */ zzq = zrefq / ( 0.620 * nu / fabs( ustar_test ) ); return zzq; /* end of subroutine z_o_zq */ } double z_o_zt( double zref, double zreft, double *ustar, double nu ) { double ustar_test, zzt; ustar_test = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] ); /* determine zreft / zT */ zzt = zreft / ( 0.40 * nu / fabs( ustar_test ) ); 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 betac_prime, double betag_prime, double betas_prime, double *betac, double *betag, double *betas, double *wind_vel, double *cp, double *wa, double *hsig, double ww, double ww_eql, int ss_prm, double ss_val, double *ustar, double *monin_inv, double *zz, int count, int i ) { short period_known; int test; double bstar, B_Toba, cmin, dum, hzzc, hzzg, period[2], p1, p2, p3, p4, phi_u, scale, sig_slope, ueff, u_hsig, Uorb, ustar_temp, wave_len, wstar2, ww_mod, z0, z0c; test = 0; ustar_temp = sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] ); cmin = sqrt( 2.0 * sqrt( G * sfcten / denwat ) ); bstar = *monin_inv * ustar_temp / KV; wstar2 = pow( fabs( sqrt( ustar[0] * ustar[0] + ustar[1] * ustar[1] ) * bstar ) * H, 0.6666667); ueff = sqrt( ( wind_vel[i] + wstar2 ) * ( wind_vel[i] + wstar2 ) ); 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( wind_vel[i] * wind_vel[i] + wstar2 ); z0 = betas[i] * 0.1100 * nu / fabs( ustar[i] ) + betac[i] * B * sfcten / ( ustar[i] * ustar[i] * denwat ) + betag[i] * 0.480 / wa[i] * ustar[i] * ustar[i] / G; phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0, monin_inv ); 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 * fabs( 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; ueff = sqrt( wind_vel[i] * wind_vel[i] + wstar2 ); z0 = betas[i] * 0.1100 * nu / fabs( ustar[i] ) + betac[i] * B * sfcten / ( ustar[i] * ustar[i] * denwat ) + betag[i] * 0.480 / wa[i] * ustar[i] * ustar[i] / G; phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0, monin_inv ); ww = ( log( hsig[i] / z0 + 1.0 ) + phi_u ) / ( KV * ww_eql * wa[i] ); /* determine the wave characteristics */ wave_len = *cp * period[i]; } else { z0 = betas[i] * 0.1100 * nu / fabs( ustar[i] ) + betac[i] * B * sfcten / ( ustar_temp * fabs( ustar[i] ) * denwat ) + betag[i] * 0.480 / wa[i] * ustar_temp * fabs( ustar[i] ) / G; /* 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 - 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; } } /* determine beta-c */ betac[i] = betac_prime; scale = 0.5; hzzg = scale * hsig[i] * zz[i] / zref; /* can't be -ve */ phi_u = phi_u_f( betag[i] * scale * hsig[i], hzzg, monin_inv ); Uorb = betag[i] * scale * 0.5 * PI2 * hsig[i] / period[i]; if ( fabs( ustar[i] ) > 0.00001 ) if ( i == 0 ) betac[i] = betac_prime * exp( -KV * Uorb / ustar[i] ); else betac[i] = betac_prime * exp( -KV * Uorb / fabs( ustar[i] ) ); else betac[i] = betac_prime; /* determine wave age */ phi_u = phi_u_f( betag[i] * hsig[i], hsig[i] / z0, monin_inv ); 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.0 */ 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] ); } else if ( !fixedwa && i == 0 ) { /* sea state determined from wave observations (other than wave age) */ if ( fabs( ustar[i] ) > 0.0 ) wa[i] = *cp / fabs( ustar[i] ); else wa[i] = 999.0; } /* 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 = wa[0] * fabs( ustar[0] ); /* determine zref / z0 */ betas[i] = betas_prime; 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.480 / wa[i] * ustar_temp * fabs( ustar[i] ) / G, 2.0 ) ) ); /* 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 ); 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 *ustar, double *qstar, double *tstar, double *monin_inv, double *zz, double *cp, double *wa, double *hsig, double CRIT, double denwat, double nu, double qmixa, double sfcten, double ta, double *wind_vel, double betac_prime, double betag_prime, double betas_prime, double *betac, double *betag, double *betas, double ww, double ww_eql, int ss_prm, double ss_val, double zref, double zreft ) { 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, qstar, tstar, 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( wind_vel[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, betac_prime, betag_prime, betas_prime, betac, betag, betas, wind_vel, cp, wa, hsig, ww, ww_eql, ss_prm, ss_val, ustar, monin_inv, zz, count, i ); phi_u = phi_u_f( zref, zz[i], monin_inv ); ustar[i] = f_ustar( fixed_uen, neutral, bstar, CON_P, phi_u, wind_vel[i], ww, zref, ustar, monin_inv, zz, 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; } }