|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! File %M% from Library %Q% 00002 ! Version %I% from %G% extracted: %H% 00003 !------------------------------------------------------------------------------ 00004 00005 MODULE mode_sfcflx 00006 00007 !------------------------------------------------------------------------------ 00008 ! 00009 ! Description: 00010 ! 00011 ! The main program unit of 00012 ! the atmospheric surface-layer parameterization scheme "sfcflx". 00013 ! "sfcflx" is used to compute fluxes 00014 ! of momentum and of sensible and latent heat over lakes. 00015 ! The surface-layer scheme developed by Mironov (1991) was used as the starting point. 00016 ! It was modified and further developed to incorporate new results as to 00017 ! the roughness lenghts for scalar quantities, 00018 ! heat and mass transfer in free convection, 00019 ! and the effect of limited fetch on the momentum transfer. 00020 ! Apart from the momentum flux and sensible and latent heat fluxes, 00021 ! the long-wave radiation flux from the water surface and 00022 ! the long-wave radiation flux from the atmosphere can also be computed. 00023 ! The atmospheric long-wave radiation flux is computed with simple empirical formulae, 00024 ! where the atmospheric emissivity is taken to be dependent on 00025 ! the water vapour pressure and cloud fraction. 00026 ! 00027 ! A description of sfcflx is available from the author. 00028 ! Dmitrii Mironov 00029 ! German Weather Service, Kaiserleistr. 29/35, D-63067 Offenbach am Main, Germany. 00030 ! dmitrii.mironov@dwd.de 00031 ! 00032 ! Lines embraced with "!_tmp" contain temporary parts of the code. 00033 ! Lines embraced/marked with "!_dev" may be replaced 00034 ! as improved parameterizations are developed and tested. 00035 ! Lines embraced/marked with "!_dm" are DM's comments 00036 ! that may be helpful to a user. 00037 ! Lines embraced/marked with "!_dbg" are used 00038 ! for debugging purposes only. 00039 ! 00040 ! 00041 ! Current Code Owner: DWD, Dmitrii Mironov 00042 ! Phone: +49-69-8062 2705 00043 ! Fax: +49-69-8062 3721 00044 ! E-mail: dmitrii.mironov@dwd.de 00045 ! 00046 ! History: 00047 ! Version Date Name 00048 ! ---------- ---------- ---- 00049 ! 1.00 2005/11/17 Dmitrii Mironov 00050 ! Initial release 00051 ! !VERSION! !DATE! <Your name> 00052 ! <Modification comments> 00053 ! 00054 ! Code Description: 00055 ! Language: Fortran 90. 00056 ! Software Standards: "European Standards for Writing and 00057 ! Documenting Exchangeable Fortran 90 Code". 00058 !============================================================================== 00059 ! 00060 ! Declarations: 00061 ! 00062 ! Modules used: 00063 00064 !USE modd_data_parameters , ONLY : & 00065 ! ireals , &! KIND-type parameter for real variables 00066 ! iintegers ! KIND-type parameter for "normal" integer variables 00067 00068 USE modd_flake_parameters , ONLY : & 00069 tpl_grav , &! Acceleration due to gravity [m s^{-2}] 00070 tpl_T_f , &! Fresh water freezing point [K] 00071 tpl_rho_w_r , &! Maximum density of fresh water [kg m^{-3}] 00072 tpl_c_w , &! Specific heat of water [J kg^{-1} K^{-1}] 00073 tpl_L_f , &! Latent heat of fusion [J kg^{-1}] 00074 h_Ice_min_flk ! Minimum ice thickness [m] 00075 00076 !============================================================================== 00077 00078 ! 00079 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00080 USE PARKIND1 ,ONLY : JPRB 00081 ! 00082 IMPLICIT NONE 00083 00084 !============================================================================== 00085 ! 00086 ! Declarations 00087 00088 ! Dimensionless constants in the Monin-Obukhov surface-layer 00089 ! similarity relations and in the expressions for the roughness lengths. 00090 REAL , PARAMETER :: 00091 c_Karman = 0.40 , ! The von Karman constant 00092 Pr_neutral = 1.0 , ! Turbulent Prandtl number at neutral static stability 00093 Sc_neutral = 1.0 , ! Turbulent Schmidt number at neutral static stability 00094 c_MO_u_stab = 5.0 , ! Constant of the MO theory (wind, stable stratification) 00095 c_MO_t_stab = 5.0 , ! Constant of the MO theory (temperature, stable stratification) 00096 c_MO_q_stab = 5.0 , ! Constant of the MO theory (humidity, stable stratification) 00097 c_MO_u_conv = 15.0 , ! Constant of the MO theory (wind, convection) 00098 c_MO_t_conv = 15.0 , ! Constant of the MO theory (temperature, convection) 00099 c_MO_q_conv = 15.0 , ! Constant of the MO theory (humidity, convection) 00100 c_MO_u_exp = 0.25 , ! Constant of the MO theory (wind, exponent) 00101 c_MO_t_exp = 0.5 , ! Constant of the MO theory (temperature, exponent) 00102 c_MO_q_exp = 0.5 , ! Constant of the MO theory (humidity, exponent) 00103 z0u_ice_rough = 1.0E-03 , ! Aerodynamic roughness of the ice surface [m] (rough flow) 00104 c_z0u_smooth = 0.1 , ! Constant in the expression for z0u (smooth flow) 00105 c_z0u_rough = 1.23E-02 , ! The Charnock constant in the expression for z0u (rough flow) 00106 c_z0u_rough_L = 1.00E-01 , ! An increased Charnock constant (used as the upper limit) 00107 c_z0u_ftch_f = 0.70 , ! Factor in the expression for fetch-dependent Charnock parameter 00108 c_z0u_ftch_ex = 0.3333333 , ! Exponent in the expression for fetch-dependent Charnock parameter 00109 c_z0t_rough_1 = 4.0 , ! Constant in the expression for z0t (factor) 00110 c_z0t_rough_2 = 3.2 , ! Constant in the expression for z0t (factor) 00111 c_z0t_rough_3 = 0.5 , ! Constant in the expression for z0t (exponent) 00112 c_z0q_rough_1 = 4.0 , ! Constant in the expression for z0q (factor) 00113 c_z0q_rough_2 = 4.2 , ! Constant in the expression for z0q (factor) 00114 c_z0q_rough_3 = 0.5 , ! Constant in the expression for z0q (exponent) 00115 c_z0t_ice_b0s = 1.250 , ! Constant in the expression for z0t over ice 00116 c_z0t_ice_b0t = 0.149 , ! Constant in the expression for z0t over ice 00117 c_z0t_ice_b1t = -0.550 , ! Constant in the expression for z0t over ice 00118 c_z0t_ice_b0r = 0.317 , ! Constant in the expression for z0t over ice 00119 c_z0t_ice_b1r = -0.565 , ! Constant in the expression for z0t over ice 00120 c_z0t_ice_b2r = -0.183 , ! Constant in the expression for z0t over ice 00121 c_z0q_ice_b0s = 1.610 , ! Constant in the expression for z0q over ice 00122 c_z0q_ice_b0t = 0.351 , ! Constant in the expression for z0q over ice 00123 c_z0q_ice_b1t = -0.628 , ! Constant in the expression for z0q over ice 00124 c_z0q_ice_b0r = 0.396 , ! Constant in the expression for z0q over ice 00125 c_z0q_ice_b1r = -0.512 , ! Constant in the expression for z0q over ice 00126 c_z0q_ice_b2r = -0.180 , ! Constant in the expression for z0q over ice 00127 Re_z0s_ice_t = 2.5 , ! Threshold value of the surface Reynolds number 00128 ! used to compute z0t and z0q over ice (Andreas 2002) 00129 Re_z0u_thresh = 0.1 ! Threshold value of the roughness Reynolds number 00130 ! [value from Zilitinkevich, Grachev, and Fairall (200), 00131 ! currently not used] 00132 00133 ! Dimensionless constants 00134 REAL , PARAMETER :: 00135 c_free_conv = 0.14 ! Constant in the expressions for fluxes in free convection 00136 00137 ! Dimensionless constants 00138 REAL , PARAMETER :: 00139 c_lwrad_emis = 0.99 ! Surface emissivity with respect to the long-wave radiation 00140 00141 ! Thermodynamic parameters 00142 REAL , PARAMETER :: 00143 tpsf_C_StefBoltz = 5.67E-08 , ! The Stefan-Boltzmann constant [W m^{-2} K^{-4}] 00144 tpsf_R_dryair = 2.8705E+02 , ! Gas constant for dry air [J kg^{-1} K^{-1}] 00145 tpsf_R_watvap = 4.6151E+02 , ! Gas constant for water vapour [J kg^{-1} K^{-1}] 00146 tpsf_c_a_p = 1.005E+03 , ! Specific heat of air at constant pressure [J kg^{-1} K^{-1}] 00147 tpsf_L_evap = 2.501E+06 , ! Specific heat of evaporation [J kg^{-1}] 00148 tpsf_nu_u_a = 1.50E-05 , ! Kinematic molecular viscosity of air [m^{2} s^{-1}] 00149 tpsf_kappa_t_a = 2.20E-05 , ! Molecular temperature conductivity of air [m^{2} s^{-1}] 00150 tpsf_kappa_q_a = 2.40E-05 ! Molecular diffusivity of air for water vapour [m^{2} s^{-1}] 00151 00152 ! Derived thermodynamic parameters 00153 REAL , PARAMETER :: 00154 tpsf_Rd_o_Rv = tpsf_R_dryair/tpsf_R_watvap , ! Ratio of gas constants (Rd/Rv) 00155 tpsf_alpha_q = (1.-tpsf_Rd_o_Rv)/tpsf_Rd_o_Rv ! Diemsnionless ratio 00156 00157 ! Thermodynamic parameters 00158 REAL , PARAMETER :: 00159 P_a_ref = 1.0E+05 ! Reference pressure [N m^{-2} = kg m^{-1} s^{-2}] 00160 00161 00162 ! The variables declared below 00163 ! are accessible to all program units of the MODULE "sfcflx" 00164 ! and to the driving routines that use "sfcflx". 00165 ! These are basically the quantities computed by sfcflx. 00166 ! Apart from these quantities, there a few local scalars 00167 ! used by sfcflx routines mainly for security reasons. 00168 ! All variables declared below have a suffix "sf". 00169 00170 ! sfcflx variables of type REAL 00171 00172 ! Roughness lengths 00173 REAL :: 00174 z0u_sf , ! Roughness length with respect to wind velocity [m] 00175 z0t_sf , ! Roughness length with respect to potential temperature [m] 00176 z0q_sf ! Roughness length with respect to specific humidity [m] 00177 00178 ! Fluxes in the surface air layer 00179 REAL :: 00180 u_star_a_sf , ! Friction velocity [m s^{-1}] 00181 Q_mom_a_sf , ! Momentum flux [N m^{-2}] 00182 Q_sens_a_sf , ! Sensible heat flux [W m^{-2}] 00183 Q_lat_a_sf , ! Laten heat flux [W m^{-2}] 00184 Q_watvap_a_sf ! Flux of water vapout [kg m^{-2} s^{-1}] 00185 00186 ! Security constants 00187 REAL , PARAMETER :: 00188 u_wind_min_sf = 1.0E-02 , ! Minimum wind speed [m s^{-1}] 00189 u_star_min_sf = 1.0E-04 , ! Minimum value of friction velocity [m s^{-1}] 00190 z0t_min_sf = 1.0E-11 , ! Minimum value of thermal roughness length [m s^{-1}] 00191 c_accur_sf = 1.0E-07 , ! A small number (accuracy) 00192 c_small_sf = 1.0E-04 ! A small number (used to compute fluxes) 00193 00194 ! Useful constants 00195 REAL , PARAMETER :: 00196 num_1o3_sf = 1./3. ! 1/3 00197 00198 !============================================================================== 00199 ! Procedures 00200 !============================================================================== 00201 00202 CONTAINS 00203 00204 !============================================================================== 00205 ! The codes of the sfcflx procedures are stored in separate "*.incf" files 00206 ! and are included below. 00207 !------------------------------------------------------------------------------ 00208 00209 !============================================================================== 00210 ! For SURFEX needs, separate *.incf files are explicitly expanded 00211 !============================================================================== 00212 !SURFEX include 'sfcflx_lwradatm.incf' 00213 ! File %M% from Library %Q% 00214 ! Version %I% from %G% extracted: %H% 00215 !------------------------------------------------------------------------------ 00216 00217 !SURFEX REAL FUNCTION sfcflx_lwradatm (T_a, e_a, cl_tot, cl_low) 00218 FUNCTION sfcflx_lwradatm (T_a, e_a, cl_tot, cl_low) 00219 00220 !------------------------------------------------------------------------------ 00221 ! 00222 ! Description: 00223 ! 00224 ! Computes the long-wave radiation flux from the atmosphere 00225 ! as function of air temperature, water vapour pressure and cloud fraction. 00226 ! 00227 ! 00228 ! Current Code Owner: DWD, Dmitrii Mironov 00229 ! Phone: +49-69-8062 2705 00230 ! Fax: +49-69-8062 3721 00231 ! E-mail: dmitrii.mironov@dwd.de 00232 ! 00233 ! History: 00234 ! Version Date Name 00235 ! ---------- ---------- ---- 00236 ! 1.00 2005/11/17 Dmitrii Mironov 00237 ! Initial release 00238 ! !VERSION! !DATE! <Your name> 00239 ! <Modification comments> 00240 ! 00241 ! Code Description: 00242 ! Language: Fortran 90. 00243 ! Software Standards: "European Standards for Writing and 00244 ! Documenting Exchangeable Fortran 90 Code". 00245 !============================================================================== 00246 ! 00247 ! Declarations: 00248 ! 00249 ! Modules used: 00250 00251 !_dm Parameters are USEd in module "sfcflx". 00252 !_nu USE modd_data_parameters , ONLY : & 00253 !_nu ireals, & ! KIND-type parameter for real variables 00254 !_nu iintegers ! KIND-type parameter for "normal" integer variables 00255 00256 !============================================================================== 00257 00258 IMPLICIT NONE 00259 00260 !============================================================================== 00261 ! 00262 ! Declarations 00263 00264 ! Input (function argument) 00265 REAL , INTENT(IN) :: 00266 T_a , ! Air temperature [K] 00267 e_a , ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 00268 cl_tot , ! Total cloud cover [0,1] 00269 cl_low ! Lowe-level cloud cover [0,1] 00270 00271 ! Output (function result) 00272 REAL :: 00273 sfcflx_lwradatm ! Long-wave radiation flux [W m^{-2}] 00274 00275 00276 ! Local parameters 00277 00278 ! Coefficients in the empirical formulation 00279 ! developed at the Main Geophysical Observatory (MGO), St. Petersburg, Russia. 00280 REAL , PARAMETER :: 00281 c_lmMGO_1 = 43.057924 , ! Empirical coefficient 00282 c_lmMGO_2 = 540.795 ! Empirical coefficient 00283 ! Temperature-dependent cloud-correction coefficients in the MGO formula 00284 INTEGER , PARAMETER :: 00285 nband_coef = 6 ! Number of temperature bands 00286 REAL , PARAMETER, DIMENSION (nband_coef) :: 00287 corr_cl_tot = (/0.70, 0.45, 0.32, 00288 0.23, 0.18, 0.13/) , ! Total clouds 00289 corr_cl_low = (/0.76, 0.49, 0.35, 00290 0.26, 0.20, 0.15/) , ! Low-level clouds 00291 corr_cl_midhigh = (/0.46, 0.30, 0.21, 00292 0.15, 0.12, 0.09/) ! Mid- and high-level clouds 00293 REAL , PARAMETER :: 00294 T_low = 253.15 , ! Low-limit temperature in the interpolation formula [K] 00295 del_T = 10.0 ! Temperature step in the interpolation formula [K] 00296 00297 ! Coefficients in the empirical water-vapour correction function 00298 ! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001). 00299 REAL , PARAMETER :: 00300 c_watvap_corr_min = 0.6100 , ! Empirical coefficient (minimum value of the correction function) 00301 c_watvap_corr_max = 0.7320 , ! Empirical coefficient (maximum value of the correction function) 00302 c_watvap_corr_e = 0.0050 ! Empirical coefficient [(N m^{-2})^{-1/2}] 00303 00304 ! Local variables of type INTEGER 00305 INTEGER :: 00306 i ! Loop index 00307 00308 ! Local variables of type REAL 00309 REAL :: 00310 c_cl_tot_corr , ! The MGO cloud correction coefficient, total clouds 00311 c_cl_low_corr , ! The MGO cloud correction coefficient, low-level clouds 00312 c_cl_midhigh_corr , ! The MGO cloud correction coefficient, mid- and high-level clouds 00313 T_corr , ! Temperature used to compute the MGO cloud correction [K] 00314 f_wvpres_corr , ! Correction function with respect to water vapour 00315 f_cloud_corr ! Correction function with respect to cloudiness 00316 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00317 00318 !============================================================================== 00319 ! Start calculations 00320 !------------------------------------------------------------------------------ 00321 00322 ! Water-vapour correction function 00323 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_LWRADATM',0,ZHOOK_HANDLE) 00324 f_wvpres_corr = c_watvap_corr_min + c_watvap_corr_e*SQRT(e_a) 00325 f_wvpres_corr = MIN(f_wvpres_corr, c_watvap_corr_max) 00326 00327 ! Cloud-correction coefficients using the MGO formulation with linear interpolation 00328 IF(T_a.LT.T_low) THEN 00329 c_cl_tot_corr = corr_cl_tot(1) 00330 c_cl_low_corr = corr_cl_low(1) 00331 c_cl_midhigh_corr = corr_cl_midhigh(1) 00332 ELSE IF(T_a.GE.T_low+(nband_coef-1)*del_T) THEN 00333 c_cl_tot_corr = corr_cl_tot(nband_coef) 00334 c_cl_low_corr = corr_cl_low(nband_coef) 00335 c_cl_midhigh_corr = corr_cl_midhigh(nband_coef) 00336 ELSE 00337 T_corr = T_low 00338 DO i=1, nband_coef-1 00339 IF(T_a.GE.T_corr.AND.T_a.LT.T_corr+del_T) THEN 00340 c_cl_tot_corr = (T_a-T_corr)/del_T 00341 c_cl_low_corr = corr_cl_low(i) + (corr_cl_low(i+1)-corr_cl_low(i))*c_cl_tot_corr 00342 c_cl_midhigh_corr = corr_cl_midhigh(i) + (corr_cl_midhigh(i+1)-corr_cl_midhigh(i))*c_cl_tot_corr 00343 c_cl_tot_corr = corr_cl_tot(i) + (corr_cl_tot(i+1)-corr_cl_tot(i))*c_cl_tot_corr 00344 END IF 00345 T_corr = T_corr + del_T 00346 END DO 00347 END IF 00348 ! Cloud correction function 00349 IF(cl_low.LT.0.) THEN ! Total cloud cover only 00350 f_cloud_corr = 1. + c_cl_tot_corr*cl_tot*cl_tot 00351 ELSE ! Total and low-level cloud cover 00352 f_cloud_corr = (1. + c_cl_low_corr*cl_low*cl_low) & 00353 * (1. + c_cl_midhigh_corr*(cl_tot*cl_tot-cl_low*cl_low)) 00354 END IF 00355 00356 ! Long-wave radiation flux [W m^{-2}] 00357 00358 ! The MGO formulation 00359 !_nu The MGO formulation 00360 !_nu sfcflx_lwradatm = -sfcflx_lwradatm*c_lwrad_emis & 00361 !_nu * (c_lmMGO_1*SQRT(tpsf_C_StefBoltz*T_a**4)-c_lmMGO_2) 00362 !_nu 00363 00364 ! "Conventional" formulation 00365 ! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001) 00366 sfcflx_lwradatm = -c_lwrad_emis*tpsf_C_StefBoltz*T_a**4 & 00367 * f_wvpres_corr*f_cloud_corr 00368 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_LWRADATM',1,ZHOOK_HANDLE) 00369 00370 !------------------------------------------------------------------------------ 00371 ! End calculations 00372 !============================================================================== 00373 00374 END FUNCTION sfcflx_lwradatm 00375 00376 00377 !============================================================================== 00378 !SURFEX include 'sfcflx_lwradwsfc.incf' 00379 ! File %M% from Library %Q% 00380 ! Version %I% from %G% extracted: %H% 00381 !------------------------------------------------------------------------------ 00382 00383 !SURFEX REAL FUNCTION sfcflx_lwradwsfc (T) 00384 FUNCTION sfcflx_lwradwsfc (emis,T) 00385 00386 !------------------------------------------------------------------------------ 00387 ! 00388 ! Description: 00389 ! 00390 ! Computes the surface long-wave radiation flux 00391 ! as function of temperature. 00392 ! 00393 ! 00394 ! Current Code Owner: DWD, Dmitrii Mironov 00395 ! Phone: +49-69-8062 2705 00396 ! Fax: +49-69-8062 3721 00397 ! E-mail: dmitrii.mironov@dwd.de 00398 ! 00399 ! History: 00400 ! Version Date Name 00401 ! ---------- ---------- ---- 00402 ! 1.00 2005/11/17 Dmitrii Mironov 00403 ! Initial release 00404 ! !VERSION! !DATE! <Your name> 00405 ! <Modification comments> 00406 ! 00407 ! Code Description: 00408 ! Language: Fortran 90. 00409 ! Software Standards: "European Standards for Writing and 00410 ! Documenting Exchangeable Fortran 90 Code". 00411 !============================================================================== 00412 ! 00413 ! Declarations: 00414 ! 00415 ! Modules used: 00416 00417 !_dm Parameters are USEd in module "sfcflx". 00418 !_nu USE modd_data_parameters , ONLY : & 00419 !_nu ireals, & ! KIND-type parameter for real variables 00420 !_nu iintegers ! KIND-type parameter for "normal" integer variables 00421 00422 !============================================================================== 00423 00424 IMPLICIT NONE 00425 00426 !============================================================================== 00427 ! 00428 ! Declarations 00429 00430 ! Input (function argument) 00431 REAL , INTENT(IN) :: 00432 emis , ! Emissivity 00433 T ! Temperature [K] 00434 00435 ! Output (function result) 00436 REAL :: 00437 sfcflx_lwradwsfc ! Long-wave radiation flux [W m^{-2}] 00438 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00439 00440 !============================================================================== 00441 ! Start calculations 00442 !------------------------------------------------------------------------------ 00443 00444 ! Long-wave radiation flux [W m^{-2}] 00445 00446 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_LWRADWSFC',0,ZHOOK_HANDLE) 00447 sfcflx_lwradwsfc = emis*tpsf_C_StefBoltz*T**4 00448 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_LWRADWSFC',1,ZHOOK_HANDLE) 00449 00450 !------------------------------------------------------------------------------ 00451 ! End calculations 00452 !============================================================================== 00453 00454 END FUNCTION sfcflx_lwradwsfc 00455 00456 00457 !============================================================================== 00458 !SURFEX include 'sfcflx_momsenlat.incf' 00459 ! File %M% from Library %Q% 00460 ! Version %I% from %G% extracted: %H% 00461 !------------------------------------------------------------------------------ 00462 00463 SUBROUTINE sfcflx_momsenlat ( height_u, height_tq, fetch, & 00464 U_a, T_a, q_a, T_s, P_a, h_ice, & 00465 Q_momentum, Q_sensible, Q_latent, Q_watvap,& 00466 Ri, z0u_ini, z0t_ini) 00467 00468 !------------------------------------------------------------------------------ 00469 ! 00470 ! Description: 00471 ! 00472 ! The sfcflx routine 00473 ! where fluxes of momentum and of sensible and latent heat 00474 ! at the air-water or air-ice (air-snow) interface are computed. 00475 ! 00476 ! Lines embraced with "!_tmp" contain temporary parts of the code. 00477 ! Lines embraced/marked with "!_dev" may be replaced 00478 ! as improved parameterizations are developed and tested. 00479 ! Lines embraced/marked with "!_dm" are DM's comments 00480 ! that may be helpful to a user. 00481 ! Lines embraced/marked with "!_dbg" are used 00482 ! for debugging purposes only. 00483 ! 00484 ! 00485 ! Current Code Owner: DWD, Dmitrii Mironov 00486 ! Phone: +49-69-8062 2705 00487 ! Fax: +49-69-8062 3721 00488 ! E-mail: dmitrii.mironov@dwd.de 00489 ! 00490 ! History: 00491 ! Version Date Name 00492 ! ---------- ---------- ---- 00493 ! 1.00 2005/11/17 Dmitrii Mironov 00494 ! Initial release 00495 ! !VERSION! !DATE! <Your name> 00496 ! <Modification comments> 00497 ! 00498 ! Code Description: 00499 ! Language: Fortran 90. 00500 ! Software Standards: "European Standards for Writing and 00501 ! Documenting Exchangeable Fortran 90 Code". 00502 !============================================================================== 00503 ! 00504 ! Declarations: 00505 ! 00506 ! Modules used: 00507 00508 !_dm Parameters are USEd in module "sfcflx". 00509 !_nu USE modd_data_parameters , ONLY : & 00510 !_nu ireals, & ! KIND-type parameter for real variables 00511 !_nu iintegers ! KIND-type parameter for "normal" integer variables 00512 00513 !============================================================================== 00514 00515 IMPLICIT NONE 00516 00517 !============================================================================== 00518 ! 00519 ! Declarations 00520 00521 ! Input (procedure arguments) 00522 00523 REAL , INTENT(IN) :: 00524 height_u , ! Height where wind is measured [m] 00525 height_tq , ! Height where temperature and humidity are measured [m] 00526 fetch , ! Typical wind fetch [m] 00527 U_a , ! Wind speed [m s^{-1}] 00528 T_a , ! Air temperature [K] 00529 q_a , ! Air specific humidity [-] 00530 T_s , ! Surface temperature (water, ice or snow) [K] 00531 P_a , ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}] 00532 h_ice , ! Ice thickness [m] 00533 ! 00534 !plm add z0u_ini and z0t_ini as input to fill z0u_sf and z0t_sf used in 00535 ! flake interface for roughness length 00536 z0u_ini , &! Initial roughness for momentum 00537 z0t_ini ! Initial roughness for heat 00538 00539 ! Output (procedure arguments) 00540 00541 REAL , INTENT(INOUT) :: 00542 Q_momentum ! Momentum flux [N m^{-2}] 00543 REAL , INTENT(OUT) :: 00544 Q_sensible , ! Sensible heat flux [W m^{-2}] 00545 Q_latent , ! Laten heat flux [W m^{-2}] 00546 Q_watvap , ! Flux of water vapout [kg m^{-2} s^{-1}] 00547 Ri ! Gradient Richardson number 00548 00549 00550 ! Local parameters of type INTEGER 00551 INTEGER , PARAMETER :: 00552 n_iter_max = 5 ! Maximum number of iterations 00553 ! n_iter_max = 24 ! Maximum number of iterations 00554 00555 ! Local variables of type LOGICAL 00556 LOGICAL :: 00557 l_conv_visc , ! Switch, TRUE = viscous free convection, the Nu=C Ra^(1/3) law is used 00558 l_conv_cbl ! Switch, TRUE = CBL scale convective structures define surface fluxes 00559 00560 ! Local variables of type INTEGER 00561 INTEGER :: 00562 i , ! Loop index 00563 n_iter ! Number of iterations performed 00564 00565 ! Local variables of type REAL 00566 REAL :: 00567 rho_a , ! Air density [kg m^{-3}] 00568 wvpres_s , ! Saturation water vapour pressure at T=T_s [N m^{-2}] 00569 q_s ! Saturation specific humidity at T=T_s [-] 00570 00571 ! Local variables of type REAL 00572 REAL :: 00573 Q_mom_tur , ! Turbulent momentum flux [N m^{-2}] 00574 Q_sen_tur , ! Turbulent sensible heat flux [W m^{-2}] 00575 Q_lat_tur , ! Turbulent laten heat flux [W m^{-2}] 00576 Q_mom_mol , ! Molecular momentum flux [N m^{-2}] 00577 Q_sen_mol , ! Molecular sensible heat flux [W m^{-2}] 00578 Q_lat_mol , ! Molecular laten heat flux [W m^{-2}] 00579 Q_mom_con , ! Momentum flux in free convection [N m^{-2}] 00580 Q_sen_con , ! Sensible heat flux in free convection [W m^{-2}] 00581 Q_lat_con ! Laten heat flux in free convection [W m^{-2}] 00582 00583 ! Local variables of type REAL 00584 REAL :: 00585 par_conv_visc , ! Viscous convection stability parameter 00586 par_conv_cbl , ! CBL convection stability parameter 00587 c_z0u_fetch , ! Fetch-dependent Charnock parameter 00588 U_a_thresh , ! Threshld value of the wind speed [m s^{-1}] 00589 u_star_thresh , ! Threshld value of friction velocity [m s^{-1}] 00590 u_star_previter , ! Friction velocity from previous iteration [m s^{-1}] 00591 u_star_n , ! Friction velocity at neutral stratification [m s^{-1}] 00592 u_star_st , ! Friction velocity with due regard for stratification [m s^{-1}] 00593 ZoL , ! The z/L ratio, z=height_u 00594 !salgado: Ri is passed as an argument (OUT) in order to be available in flake_interface 00595 !Ri , & ! Gradient Richardson number 00596 Ri_cr , &! Critical value of Ri 00597 R_z , &! Ratio of "height_tq" to "height_u" 00598 Fun , &! A function of generic variable "x" 00599 Fun_prime , &! Derivative of "Fun" with respect to "x" 00600 Delta , &! Relative error 00601 psi_u , &! The MO stability function for wind profile 00602 psi_t , &! The MO stability function for temperature profile 00603 psi_q ! The MO stability function for specific humidity profile 00604 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00605 00606 00607 !============================================================================== 00608 ! Start calculations 00609 !------------------------------------------------------------------------------ 00610 00611 !_dm All fluxes are positive when directed upwards. 00612 00613 !------------------------------------------------------------------------------ 00614 ! Compute saturation specific humidity and the air density at T=T_s 00615 !------------------------------------------------------------------------------ 00616 00617 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_MOMSENLAT',0,ZHOOK_HANDLE) 00618 wvpres_s = sfcflx_satwvpres(T_s, h_ice) ! Saturation water vapour pressure at T=T_s 00619 q_s = sfcflx_spechum (wvpres_s, P_a) ! Saturation specific humidity at T=T_s 00620 rho_a = sfcflx_rhoair(T_s, q_s, P_a) ! Air density at T_s and q_s (surface values) 00621 00622 !------------------------------------------------------------------------------ 00623 ! Compute molecular fluxes of momentum and of sensible and latent heat 00624 !------------------------------------------------------------------------------ 00625 00626 !_dm The fluxes are in kinematic units 00627 Q_mom_mol = -tpsf_nu_u_a*U_a/height_u 00628 Q_sen_mol = -tpsf_kappa_t_a*(T_a-T_s)/height_tq 00629 Q_lat_mol = -tpsf_kappa_q_a*(q_a-q_s)/height_tq 00630 00631 !------------------------------------------------------------------------------ 00632 ! Compute fluxes in free convection 00633 !------------------------------------------------------------------------------ 00634 00635 par_conv_visc = (T_s-T_a)/T_s*SQRT(tpsf_kappa_t_a) + (q_s-q_a)*tpsf_alpha_q*SQRT(tpsf_kappa_q_a) 00636 IF(par_conv_visc.GT.0.) THEN ! Viscous convection takes place 00637 l_conv_visc = .TRUE. 00638 par_conv_visc = (par_conv_visc*tpl_grav/tpsf_nu_u_a)**num_1o3_sf 00639 Q_sen_con = c_free_conv*SQRT(tpsf_kappa_t_a)*par_conv_visc 00640 Q_sen_con = Q_sen_con*(T_s-T_a) 00641 Q_lat_con = c_free_conv*SQRT(tpsf_kappa_q_a)*par_conv_visc 00642 Q_lat_con = Q_lat_con*(q_s-q_a) 00643 ELSE ! No viscous convection, set fluxes to zero 00644 l_conv_visc = .FALSE. 00645 Q_sen_con = 0. 00646 Q_lat_con = 0. 00647 END IF 00648 Q_mom_con = 0. ! Momentum flux in free (viscous or CBL-scale) convection is zero 00649 00650 !------------------------------------------------------------------------------ 00651 ! Compute turbulent fluxes 00652 !------------------------------------------------------------------------------ 00653 00654 R_z = height_tq/height_u ! Ratio of "height_tq" to "height_u" 00655 Ri_cr = c_MO_t_stab/c_MO_u_stab**2*R_z ! Critical Ri 00656 Ri = tpl_grav*((T_a-T_s)/T_s+tpsf_alpha_q*(q_a-q_s))/MAX(U_a,u_wind_min_sf)**2 00657 Ri = Ri*height_u/Pr_neutral ! Gradient Richardson number 00658 00659 Turb_Fluxes: IF(U_a.LT.u_wind_min_sf.OR.Ri.GT.Ri_cr-c_small_sf) THEN ! Low wind or Ri>Ri_cr 00660 00661 u_star_st = 0. ! Set turbulent fluxes to zero 00662 Q_mom_tur = 0. 00663 Q_sen_tur = 0. 00664 Q_lat_tur = 0. 00665 z0u_sf = z0u_ini 00666 z0t_sf = z0t_ini 00667 00668 ELSE Turb_Fluxes ! Compute turbulent fluxes using MO similarity 00669 00670 ! Compute z/L, where z=height_u 00671 IF(Ri.GE.0.) THEN ! Stable stratification 00672 ZoL = SQRT(1.-4.*(c_MO_u_stab-R_z*c_MO_t_stab)*Ri) 00673 ZoL = ZoL - 1. + 2.*c_MO_u_stab*Ri 00674 ZoL = ZoL/2./c_MO_u_stab/c_MO_u_stab/(Ri_cr-Ri) 00675 ELSE ! Convection 00676 n_iter = 0 00677 Delta = 1. ! Set initial error to a large value (as compared to the accuracy) 00678 u_star_previter = Ri*MAX(1., SQRT(R_z*c_MO_t_conv/c_MO_u_conv)) ! Initial guess for ZoL 00679 DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max) 00680 Fun = u_star_previter**2*(c_MO_u_conv*u_star_previter-1.) & 00681 + Ri**2*(1.-R_z*c_MO_t_conv*u_star_previter) 00682 Fun_prime = 3.*c_MO_u_conv*u_star_previter**2 & 00683 - 2.*u_star_previter - R_z*c_MO_t_conv*Ri**2 00684 ZoL = u_star_previter - Fun/Fun_prime 00685 Delta = ABS(ZoL-u_star_previter)/MAX(c_accur_sf, ABS(ZoL+u_star_previter)) 00686 u_star_previter = ZoL 00687 n_iter = n_iter + 1 00688 END DO 00689 !_dbg 00690 ! IF(n_iter.GE.n_iter_max-1) & 00691 ! WRITE(*,*) 'ZoL: Max No. iters. exceeded (n_iter = ', n_iter, ')!' 00692 !_dbg 00693 END IF 00694 00695 ! Compute fetch-dependent Charnock parameter, use "u_star_min_sf" 00696 CALL sfcflx_roughness (fetch, U_a, u_star_min_sf, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf) 00697 00698 ! Threshold value of wind speed 00699 u_star_st = u_star_thresh 00700 CALL sfcflx_roughness (fetch, U_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf) 00701 IF(ZoL.GT.0.) THEN ! MO function in stable stratification 00702 psi_u = c_MO_u_stab*ZoL*(1.-MIN(z0u_sf/height_u, 1.)) 00703 ELSE ! MO function in convection 00704 psi_t = (1.-c_MO_u_conv*ZoL)**c_MO_u_exp 00705 psi_q = (1.-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.))**c_MO_u_exp 00706 psi_u = 2.*(ATAN(psi_t)-ATAN(psi_q)) & 00707 + 2.*LOG((1.+psi_q)/(1.+psi_t)) & 00708 + LOG((1.+psi_q*psi_q)/(1.+psi_t*psi_t)) 00709 END IF 00710 U_a_thresh = u_star_thresh/c_Karman*(LOG(height_u/z0u_sf)+psi_u) 00711 00712 ! Compute friction velocity 00713 n_iter = 0 00714 Delta = 1. ! Set initial error to a large value (as compared to the accuracy) 00715 !u_star_previter = u_star_thresh ! Initial guess for friction velocity 00716 !* modif. V. Masson (Meteo-France) : uses previous time-step momentum flux for ustar guess 00717 u_star_previter = max ( sqrt( - Q_momentum / rho_a ) , u_star_thresh ) 00718 ! 00719 IF(U_a.LE.U_a_thresh) THEN ! Smooth surface 00720 DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max) 00721 CALL sfcflx_roughness (fetch, U_a, MIN(u_star_thresh, u_star_previter), h_ice, & 00722 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf) 00723 IF(ZoL.GE.0.) THEN ! Stable stratification 00724 psi_u = c_MO_u_stab*ZoL*(1.-MIN(z0u_sf/height_u, 1.)) 00725 Fun = LOG(height_u/z0u_sf) + psi_u 00726 Fun_prime = (Fun + 1. + c_MO_u_stab*ZoL*MIN(z0u_sf/height_u, 1.))/c_Karman 00727 Fun = Fun*u_star_previter/c_Karman - U_a 00728 ELSE ! Convection 00729 psi_t = (1.-c_MO_u_conv*ZoL)**c_MO_u_exp 00730 psi_q = (1.-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.))**c_MO_u_exp 00731 psi_u = 2.*(ATAN(psi_t)-ATAN(psi_q)) & 00732 + 2.*LOG((1.+psi_q)/(1.+psi_t)) & 00733 + LOG((1.+psi_q*psi_q)/(1.+psi_t*psi_t)) 00734 Fun = LOG(height_u/z0u_sf) + psi_u 00735 Fun_prime = (Fun + 1./psi_q)/c_Karman 00736 Fun = Fun*u_star_previter/c_Karman - U_a 00737 END IF 00738 u_star_st = u_star_previter - Fun/Fun_prime 00739 Delta = ABS((u_star_st-u_star_previter)/(u_star_st+u_star_previter)) 00740 u_star_previter = u_star_st 00741 n_iter = n_iter + 1 00742 END DO 00743 ELSE ! Rough surface 00744 DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max.AND.z0t_sf>z0t_min_sf) 00745 CALL sfcflx_roughness (fetch, U_a, MAX(u_star_thresh, u_star_previter), h_ice, & 00746 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf) 00747 IF(ZoL.GE.0.) THEN ! Stable stratification 00748 psi_u = c_MO_u_stab*ZoL*(1.-MIN(z0u_sf/height_u, 1.)) 00749 Fun = LOG(height_u/z0u_sf) + psi_u 00750 Fun_prime = (Fun - 2. - 2.*c_MO_u_stab*ZoL*MIN(z0u_sf/height_u, 1.))/c_Karman 00751 Fun = Fun*u_star_previter/c_Karman - U_a 00752 ELSE ! Convection 00753 psi_t = (1.-c_MO_u_conv*ZoL)**c_MO_u_exp 00754 psi_q = (1.-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.))**c_MO_u_exp 00755 psi_u = 2.*(ATAN(psi_t)-ATAN(psi_q)) & 00756 + 2.*LOG((1.+psi_q)/(1.+psi_t)) & 00757 + LOG((1.+psi_q*psi_q)/(1.+psi_t*psi_t)) 00758 Fun = LOG(height_u/z0u_sf) + psi_u 00759 Fun_prime = (Fun - 2./psi_q)/c_Karman 00760 Fun = Fun*u_star_previter/c_Karman - U_a 00761 END IF 00762 IF(h_ice.GE.h_Ice_min_flk) THEN ! No iteration is required for rough flow over ice 00763 u_star_st = c_Karman*U_a/MAX(c_small_sf, LOG(height_u/z0u_sf)+psi_u) 00764 u_star_previter = u_star_st 00765 ELSE ! Iterate in case of open water 00766 u_star_st = u_star_previter - Fun/Fun_prime 00767 END IF 00768 Delta = ABS((u_star_st-u_star_previter)/(u_star_st+u_star_previter)) 00769 u_star_previter = u_star_st 00770 n_iter = n_iter + 1 00771 END DO 00772 END IF 00773 ! 00774 !* Modification: V. Masson (Meteo-France) 00775 !* in case convergence did not occur : 00776 ! - one keeps the values of the previous time-step for momentum roughness length 00777 ! - one recomputes the thermal and water vapor roughness lengthes 00778 ! - one estimates the momentum flux using neutral law 00779 IF (n_iter == n_iter_max .OR. z0T_sf<=z0t_min_sf) THEN 00780 u_star_st = max ( sqrt( - Q_momentum / rho_a ) , u_star_min_sf ) 00781 CALL sfcflx_roughness (fetch, U_a, u_star_st, h_ice, & 00782 c_z0u_fetch, u_star_previter, z0u_sf, z0t_sf, z0q_sf) 00783 z0u_sf = z0u_ini 00784 u_star_st = c_Karman*U_a/MAX(c_small_sf, LOG(height_u/z0u_sf)) 00785 END IF 00786 00787 !_dbg 00788 ! WRITE(*,*) 'MO stab. func. psi_u = ', psi_u, ' n_iter = ', n_iter 00789 ! WRITE(*,*) ' Wind speed = ', U_a, ' u_* = ', u_star_st 00790 ! WRITE(*,*) ' Fun = ', Fun 00791 !_dbg 00792 00793 !_dbg 00794 ! IF(n_iter.GE.n_iter_max-1) & 00795 ! WRITE(*,*) 'u_*: Max No. iters. exceeded (n_iter = ', n_iter, ')!' 00796 !_dbg 00797 00798 ! Momentum flux 00799 Q_mom_tur = -u_star_st*u_star_st 00800 00801 ! Temperature and specific humidity fluxes 00802 CALL sfcflx_roughness (fetch, U_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf) 00803 IF(ZoL.GE.0.) THEN ! Stable stratification 00804 psi_t = c_MO_t_stab*R_z*ZoL*(1.-MIN(z0t_sf/height_tq, 1.)) 00805 psi_q = c_MO_q_stab*R_z*ZoL*(1.-MIN(z0q_sf/height_tq, 1.)) 00806 !_dbg 00807 ! WRITE(*,*) 'STAB: psi_t = ', psi_t, ' psi_q = ', psi_q 00808 !_dbg 00809 ELSE ! Convection 00810 psi_u = (1.-c_MO_t_conv*R_z*ZoL)**c_MO_t_exp 00811 psi_t = (1.-c_MO_t_conv*R_z*ZoL*MIN(z0t_sf/height_tq, 1.))**c_MO_t_exp 00812 psi_t = 2.*LOG((1.+psi_t)/(1.+psi_u)) 00813 psi_u = (1.-c_MO_q_conv*R_z*ZoL)**c_MO_q_exp 00814 psi_q = (1.-c_MO_q_conv*R_z*ZoL*MIN(z0q_sf/height_tq, 1.))**c_MO_q_exp 00815 psi_q = 2.*LOG((1.+psi_q)/(1.+psi_u)) 00816 !_dbg 00817 ! WRITE(*,*) 'CONV: psi_t = ', psi_t, ' psi_q = ', psi_q 00818 !_dbg 00819 END IF 00820 Q_sen_tur = -(T_a-T_s)*u_star_st*c_Karman/Pr_neutral & 00821 / MAX(c_small_sf, LOG(height_tq/z0t_sf)+psi_t) 00822 Q_lat_tur = -(q_a-q_s)*u_star_st*c_Karman/Sc_neutral & 00823 / MAX(c_small_sf, LOG(height_tq/z0q_sf)+psi_q) 00824 00825 END IF Turb_Fluxes 00826 00827 !------------------------------------------------------------------------------ 00828 ! Decide between turbulent, molecular, and convective fluxes 00829 !------------------------------------------------------------------------------ 00830 00831 Q_momentum = MIN(Q_mom_tur, Q_mom_mol, Q_mom_con) ! Momentum flux is negative 00832 IF(l_conv_visc) THEN ! Convection, take fluxes that are maximal in magnitude 00833 IF(ABS(Q_sen_tur).GE.ABS(Q_sen_con)) THEN 00834 Q_sensible = Q_sen_tur 00835 ELSE 00836 Q_sensible = Q_sen_con 00837 END IF 00838 IF(ABS(Q_sensible).LT.ABS(Q_sen_mol)) THEN 00839 Q_sensible = Q_sen_mol 00840 END IF 00841 IF(ABS(Q_lat_tur).GE.ABS(Q_lat_con)) THEN 00842 Q_latent = Q_lat_tur 00843 ELSE 00844 Q_latent = Q_lat_con 00845 END IF 00846 IF(ABS(Q_latent).LT.ABS(Q_lat_mol)) THEN 00847 Q_latent = Q_lat_mol 00848 END IF 00849 ELSE ! Stable or neutral stratification, chose fluxes that are maximal in magnitude 00850 IF(ABS(Q_sen_tur).GE.ABS(Q_sen_mol)) THEN 00851 Q_sensible = Q_sen_tur 00852 ELSE 00853 Q_sensible = Q_sen_mol 00854 END IF 00855 IF(ABS(Q_lat_tur).GE.ABS(Q_lat_mol)) THEN 00856 Q_latent = Q_lat_tur 00857 ELSE 00858 Q_latent = Q_lat_mol 00859 END IF 00860 END IF 00861 00862 !------------------------------------------------------------------------------ 00863 ! Set output (notice that fluxes are no longer in kinematic units) 00864 !------------------------------------------------------------------------------ 00865 00866 Q_momentum = Q_momentum*rho_a 00867 Q_sensible = Q_sensible*rho_a*tpsf_c_a_p 00868 Q_watvap = Q_latent*rho_a 00869 Q_latent = tpsf_L_evap 00870 IF(h_ice.GE.h_Ice_min_flk) Q_latent = Q_latent + tpl_L_f ! Add latent heat of fusion over ice 00871 Q_latent = Q_watvap*Q_latent 00872 00873 ! Set "*_sf" variables to make fluxes accessible to driving routines that use "sfcflx" 00874 u_star_a_sf = u_star_st 00875 Q_mom_a_sf = Q_momentum 00876 Q_sens_a_sf = Q_sensible 00877 Q_lat_a_sf = Q_latent 00878 Q_watvap_a_sf = Q_watvap 00879 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_MOMSENLAT',1,ZHOOK_HANDLE) 00880 00881 !------------------------------------------------------------------------------ 00882 ! End calculations 00883 !============================================================================== 00884 00885 END SUBROUTINE sfcflx_momsenlat 00886 00887 00888 !============================================================================== 00889 00890 !============================================================================== 00891 !SURFEX include 'sfcflx_rhoair.incf' 00892 ! File %M% from Library %Q% 00893 ! Version %I% from %G% extracted: %H% 00894 !------------------------------------------------------------------------------ 00895 00896 !SURFEX REAL FUNCTION sfcflx_rhoair (T, q, P) 00897 FUNCTION sfcflx_rhoair (T, q, P) 00898 00899 !------------------------------------------------------------------------------ 00900 ! 00901 ! Description: 00902 ! 00903 ! Computes the air density as function 00904 ! of temperature, specific humidity and pressure. 00905 ! 00906 ! 00907 ! Current Code Owner: DWD, Dmitrii Mironov 00908 ! Phone: +49-69-8062 2705 00909 ! Fax: +49-69-8062 3721 00910 ! E-mail: dmitrii.mironov@dwd.de 00911 ! 00912 ! History: 00913 ! Version Date Name 00914 ! ---------- ---------- ---- 00915 ! 1.00 2005/11/17 Dmitrii Mironov 00916 ! Initial release 00917 ! !VERSION! !DATE! <Your name> 00918 ! <Modification comments> 00919 ! 00920 ! Code Description: 00921 ! Language: Fortran 90. 00922 ! Software Standards: "European Standards for Writing and 00923 ! Documenting Exchangeable Fortran 90 Code". 00924 !============================================================================== 00925 ! 00926 ! Declarations: 00927 ! 00928 ! Modules used: 00929 00930 !_dm Parameters are USEd in module "sfcflx". 00931 !_nu USE modd_data_parameters , ONLY : & 00932 !_nu ireals, & ! KIND-type parameter for real variables 00933 !_nu iintegers ! KIND-type parameter for "normal" integer variables 00934 00935 !============================================================================== 00936 00937 IMPLICIT NONE 00938 00939 !============================================================================== 00940 ! 00941 ! Declarations 00942 00943 ! Input (function argument) 00944 REAL , INTENT(IN) :: 00945 T , ! Temperature [K] 00946 q , ! Specific humidity 00947 P ! Pressure [N m^{-2} = kg m^{-1} s^{-2}] 00948 00949 ! Output (function result) 00950 REAL :: 00951 sfcflx_rhoair ! Air density [kg m^{-3}] 00952 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00953 00954 !============================================================================== 00955 ! Start calculations 00956 !------------------------------------------------------------------------------ 00957 00958 ! Air density [kg m^{-3}] 00959 00960 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_RHOAIR',0,ZHOOK_HANDLE) 00961 sfcflx_rhoair = P/tpsf_R_dryair/T/(1.+(1./tpsf_Rd_o_Rv-1.)*q) 00962 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_RHOAIR',1,ZHOOK_HANDLE) 00963 00964 !------------------------------------------------------------------------------ 00965 ! End calculations 00966 !============================================================================== 00967 00968 END FUNCTION sfcflx_rhoair 00969 00970 00971 !============================================================================== 00972 !SURFEX include 'sfcflx_roughness.incf' 00973 ! File %M% from Library %Q% 00974 ! Version %I% from %G% extracted: %H% 00975 !------------------------------------------------------------------------------ 00976 00977 SUBROUTINE sfcflx_roughness (fetch, U_a, u_star, h_ice, & 00978 c_z0u_fetch, u_star_thresh, z0u, z0t, z0q) 00979 00980 !------------------------------------------------------------------------------ 00981 ! 00982 ! Description: 00983 ! 00984 ! Computes the water-surface or the ice-surface roughness lengths 00985 ! with respect to wind velocity, potential temperature and specific humidity. 00986 ! 00987 ! The water-surface roughness lengths with respect to wind velocity is computed 00988 ! from the Charnock formula when the surface is aerodynamically rough. 00989 ! A simple empirical formulation is used to account for the dependence 00990 ! of the Charnock parameter on the wind fetch. 00991 ! When the flow is aerodynamically smooth, the roughness length with respect to 00992 ! wind velocity is proportional to the depth of the viscous sub-layer. 00993 ! The water-surface roughness lengths for scalars are computed using the power-law 00994 ! formulations in terms of the roughness Reynolds number (Zilitinkevich et al. 2001). 00995 ! The ice-surface aerodynamic roughness is taken to be constant. 00996 ! The ice-surface roughness lengths for scalars 00997 ! are computed through the power-law formulations 00998 ! in terms of the roughness Reynolds number (Andreas 2002). 00999 ! 01000 ! 01001 ! Current Code Owner: DWD, Dmitrii Mironov 01002 ! Phone: +49-69-8062 2705 01003 ! Fax: +49-69-8062 3721 01004 ! E-mail: dmitrii.mironov@dwd.de 01005 ! 01006 ! History: 01007 ! Version Date Name 01008 ! ---------- ---------- ---- 01009 ! 1.00 2005/11/17 Dmitrii Mironov 01010 ! Initial release 01011 ! !VERSION! !DATE! <Your name> 01012 ! <Modification comments> 01013 ! 01014 ! Code Description: 01015 ! Language: Fortran 90. 01016 ! Software Standards: "European Standards for Writing and 01017 ! Documenting Exchangeable Fortran 90 Code". 01018 !============================================================================== 01019 ! 01020 ! Declarations: 01021 ! 01022 ! Modules used: 01023 01024 !_dm Parameters are USEd in module "sfcflx". 01025 !_nu USE modd_data_parameters , ONLY : & 01026 !_nu ireals , & ! KIND-type parameter for real variables 01027 !_nu iintegers ! KIND-type parameter for "normal" integer variables 01028 01029 !============================================================================== 01030 01031 IMPLICIT NONE 01032 01033 !============================================================================== 01034 ! 01035 ! Declarations 01036 01037 ! Input (procedure arguments) 01038 REAL , INTENT(IN) :: 01039 fetch , ! Typical wind fetch [m] 01040 U_a , ! Wind speed [m s^{-1}] 01041 u_star , ! Friction velocity in the surface air layer [m s^{-1}] 01042 h_ice ! Ice thickness [m] 01043 01044 ! Output (procedure arguments) 01045 REAL , INTENT(OUT) :: 01046 c_z0u_fetch , ! Fetch-dependent Charnock parameter 01047 u_star_thresh , ! Threshold value of friction velocity [m s^{-1}] 01048 z0u , ! Roughness length with respect to wind velocity [m] 01049 z0t , ! Roughness length with respect to potential temperature [m] 01050 z0q ! Roughness length with respect to specific humidity [m] 01051 01052 ! Local variables of type REAL 01053 REAL :: 01054 Re_s , ! Surface Reynolds number 01055 Re_s_thresh ! Threshold value of Re_s 01056 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01057 01058 !============================================================================== 01059 ! Start calculations 01060 !------------------------------------------------------------------------------ 01061 01062 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_ROUGHNESS',0,ZHOOK_HANDLE) 01063 Water_or_Ice: IF(h_ice.LT.h_Ice_min_flk) THEN ! Water surface 01064 01065 ! The Charnock parameter as dependent on dimensionless fetch 01066 c_z0u_fetch = MAX(U_a, u_wind_min_sf)**2/tpl_grav/fetch ! Inverse dimensionless fetch 01067 c_z0u_fetch = c_z0u_rough + c_z0u_ftch_f*c_z0u_fetch**c_z0u_ftch_ex 01068 c_z0u_fetch = MIN(c_z0u_fetch, c_z0u_rough_L) ! Limit Charnock parameter 01069 01070 ! Threshold value of friction velocity 01071 u_star_thresh = (c_z0u_smooth/c_z0u_fetch*tpl_grav*tpsf_nu_u_a)**num_1o3_sf 01072 01073 ! Surface Reynolds number and its threshold value 01074 Re_s = u_star**3/tpsf_nu_u_a/tpl_grav 01075 Re_s_thresh = c_z0u_smooth/c_z0u_fetch 01076 01077 ! Aerodynamic roughness 01078 IF(Re_s.LE.Re_s_thresh) THEN 01079 z0u = c_z0u_smooth*tpsf_nu_u_a/u_star ! Smooth flow 01080 ELSE 01081 z0u = c_z0u_fetch*u_star*u_star/tpl_grav ! Rough flow 01082 END IF 01083 ! Roughness for scalars 01084 z0q = c_z0u_fetch*MAX(Re_s, Re_s_thresh) 01085 z0t = c_z0t_rough_1*z0q**c_z0t_rough_3 - c_z0t_rough_2 01086 z0q = c_z0q_rough_1*z0q**c_z0q_rough_3 - c_z0q_rough_2 01087 z0t = z0u*EXP(-c_Karman/Pr_neutral*z0t) 01088 z0q = z0u*EXP(-c_Karman/Sc_neutral*z0q) 01089 01090 ELSE Water_or_Ice ! Ice surface 01091 01092 ! The Charnock parameter is not used over ice, formally set "c_z0u_fetch" to its minimum value 01093 c_z0u_fetch = c_z0u_rough 01094 01095 ! Threshold value of friction velocity 01096 u_star_thresh = c_z0u_smooth*tpsf_nu_u_a/z0u_ice_rough 01097 01098 ! Aerodynamic roughness 01099 z0u = MAX(z0u_ice_rough, c_z0u_smooth*tpsf_nu_u_a/u_star) 01100 01101 ! Roughness Reynolds number 01102 Re_s = MAX(u_star*z0u/tpsf_nu_u_a, c_accur_sf) 01103 01104 ! Roughness for scalars 01105 IF(Re_s.LE.Re_z0s_ice_t) THEN 01106 z0t = c_z0t_ice_b0t + c_z0t_ice_b1t*LOG(Re_s) 01107 z0t = MIN(z0t, c_z0t_ice_b0s) 01108 z0q = c_z0q_ice_b0t + c_z0q_ice_b1t*LOG(Re_s) 01109 z0q = MIN(z0q, c_z0q_ice_b0s) 01110 ELSE 01111 z0t = c_z0t_ice_b0r + c_z0t_ice_b1r*LOG(Re_s) + c_z0t_ice_b2r*LOG(Re_s)**2 01112 z0q = c_z0q_ice_b0r + c_z0q_ice_b1r*LOG(Re_s) + c_z0q_ice_b2r*LOG(Re_s)**2 01113 END IF 01114 z0t = z0u*EXP(z0t) 01115 z0q = z0u*EXP(z0q) 01116 01117 END IF Water_or_Ice 01118 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_ROUGHNESS',1,ZHOOK_HANDLE) 01119 01120 !------------------------------------------------------------------------------ 01121 ! End calculations 01122 !============================================================================== 01123 01124 END SUBROUTINE sfcflx_roughness 01125 01126 !============================================================================== 01127 !SURFEX include 'sfcflx_satwvpres.incf' 01128 ! File %M% from Library %Q% 01129 ! Version %I% from %G% extracted: %H% 01130 !------------------------------------------------------------------------------ 01131 01132 !SURFEX REAL FUNCTION sfcflx_satwvpres (T, h_ice) 01133 FUNCTION sfcflx_satwvpres (T, h_ice) 01134 01135 !------------------------------------------------------------------------------ 01136 ! 01137 ! Description: 01138 ! 01139 ! Computes saturation water vapour pressure 01140 ! over the water surface or over the ice surface 01141 ! as function of temperature. 01142 ! 01143 ! 01144 ! Current Code Owner: DWD, Dmitrii Mironov 01145 ! Phone: +49-69-8062 2705 01146 ! Fax: +49-69-8062 3721 01147 ! E-mail: dmitrii.mironov@dwd.de 01148 ! 01149 ! History: 01150 ! Version Date Name 01151 ! ---------- ---------- ---- 01152 ! 1.00 2005/11/17 Dmitrii Mironov 01153 ! Initial release 01154 ! !VERSION! !DATE! <Your name> 01155 ! <Modification comments> 01156 ! 01157 ! Code Description: 01158 ! Language: Fortran 90. 01159 ! Software Standards: "European Standards for Writing and 01160 ! Documenting Exchangeable Fortran 90 Code". 01161 !============================================================================== 01162 ! 01163 ! Declarations: 01164 ! 01165 ! Modules used: 01166 01167 !_dm Parameters are USEd in module "sfcflx". 01168 !_nu USE modd_data_parameters , ONLY : & 01169 !_nu ireals, & ! KIND-type parameter for real variables 01170 !_nu iintegers ! KIND-type parameter for "normal" integer variables 01171 01172 !_dm The variable is USEd in module "sfcflx". 01173 !_nu USE flake_parameters , ONLY : & 01174 !_nu h_Ice_min_flk ! Minimum ice thickness [m] 01175 01176 !============================================================================== 01177 01178 IMPLICIT NONE 01179 01180 !============================================================================== 01181 ! 01182 ! Declarations 01183 01184 ! Input (function argument) 01185 REAL , INTENT(IN) :: 01186 T , ! Temperature [K] 01187 h_ice ! Ice thickness [m] 01188 01189 ! Output (function result) 01190 REAL :: 01191 sfcflx_satwvpres ! Saturation water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 01192 01193 ! Local parameters 01194 REAL , PARAMETER :: 01195 b1_vap = 610.78 , ! Coefficient [N m^{-2} = kg m^{-1} s^{-2}] 01196 b3_vap = 273.16 , ! Triple point [K] 01197 b2w_vap = 17.2693882 , ! Coefficient (water) 01198 b2i_vap = 21.8745584 , ! Coefficient (ice) 01199 b4w_vap = 35.86 , ! Coefficient (temperature) [K] 01200 b4i_vap = 7.66 ! Coefficient (temperature) [K] 01201 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01202 01203 !============================================================================== 01204 ! Start calculations 01205 !------------------------------------------------------------------------------ 01206 01207 ! Saturation water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 01208 01209 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_SATWVPRES',0,ZHOOK_HANDLE) 01210 IF(h_ice.LT.h_Ice_min_flk) THEN ! Water surface 01211 sfcflx_satwvpres = b1_vap*EXP(b2w_vap*(T-b3_vap)/(T-b4w_vap)) 01212 ELSE ! Ice surface 01213 sfcflx_satwvpres = b1_vap*EXP(b2i_vap*(T-b3_vap)/(T-b4i_vap)) 01214 END IF 01215 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_SATWVPRES',1,ZHOOK_HANDLE) 01216 01217 !------------------------------------------------------------------------------ 01218 ! End calculations 01219 !============================================================================== 01220 01221 END FUNCTION sfcflx_satwvpres 01222 01223 01224 !============================================================================== 01225 !SURFEX include 'sfcflx_spechum.incf' 01226 ! File %M% from Library %Q% 01227 ! Version %I% from %G% extracted: %H% 01228 !------------------------------------------------------------------------------ 01229 01230 !SURFEX REAL FUNCTION sfcflx_spechum (wvpres, P) 01231 FUNCTION sfcflx_spechum (wvpres, P) 01232 01233 !------------------------------------------------------------------------------ 01234 ! 01235 ! Description: 01236 ! 01237 ! Computes specific humidity as function 01238 ! of water vapour pressure and air pressure. 01239 ! 01240 ! 01241 ! Current Code Owner: DWD, Dmitrii Mironov 01242 ! Phone: +49-69-8062 2705 01243 ! Fax: +49-69-8062 3721 01244 ! E-mail: dmitrii.mironov@dwd.de 01245 ! 01246 ! History: 01247 ! Version Date Name 01248 ! ---------- ---------- ---- 01249 ! 1.00 2005/11/17 Dmitrii Mironov 01250 ! Initial release 01251 ! !VERSION! !DATE! <Your name> 01252 ! <Modification comments> 01253 ! 01254 ! Code Description: 01255 ! Language: Fortran 90. 01256 ! Software Standards: "European Standards for Writing and 01257 ! Documenting Exchangeable Fortran 90 Code". 01258 !============================================================================== 01259 ! 01260 ! Declarations: 01261 ! 01262 ! Modules used: 01263 01264 !_dm Parameters are USEd in module "sfcflx". 01265 !_nu USE modd_data_parameters , ONLY : & 01266 !_nu ireals, & ! KIND-type parameter for real variables 01267 !_nu iintegers ! KIND-type parameter for "normal" integer variables 01268 01269 !============================================================================== 01270 01271 IMPLICIT NONE 01272 01273 !============================================================================== 01274 ! 01275 ! Declarations 01276 01277 ! Input (function argument) 01278 REAL , INTENT(IN) :: 01279 wvpres , ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 01280 P ! Air pressure [N m^{-2} = kg m^{-1} s^{-2}] 01281 01282 ! Output (function result) 01283 REAL :: 01284 sfcflx_spechum ! Specific humidity 01285 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01286 01287 !============================================================================== 01288 ! Start calculations 01289 !------------------------------------------------------------------------------ 01290 01291 ! Specific humidity 01292 01293 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_SPECHUM',0,ZHOOK_HANDLE) 01294 sfcflx_spechum = tpsf_Rd_o_Rv*wvpres/(P-(1.-tpsf_Rd_o_Rv)*wvpres) 01295 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_SPECHUM',1,ZHOOK_HANDLE) 01296 01297 !------------------------------------------------------------------------------ 01298 ! End calculations 01299 !============================================================================== 01300 01301 END FUNCTION sfcflx_spechum 01302 01303 01304 !============================================================================== 01305 !SURFEX include 'sfcflx_wvpreswetbulb.incf' 01306 ! File %M% from Library %Q% 01307 ! Version %I% from %G% extracted: %H% 01308 !------------------------------------------------------------------------------ 01309 01310 !SURFEX REAL FUNCTION sfcflx_wvpreswetbulb (T_dry, T_wetbulb, satwvpres_bulb, P) 01311 FUNCTION sfcflx_wvpreswetbulb (T_dry, T_wetbulb, satwvpres_bulb, P) 01312 01313 !------------------------------------------------------------------------------ 01314 ! 01315 ! Description: 01316 ! 01317 ! Computes water vapour pressure as function of air temperature, 01318 ! wet bulb temperature, satururation vapour pressure at wet-bulb temperature, 01319 ! and air pressure. 01320 ! 01321 ! 01322 ! Current Code Owner: DWD, Dmitrii Mironov 01323 ! Phone: +49-69-8062 2705 01324 ! Fax: +49-69-8062 3721 01325 ! E-mail: dmitrii.mironov@dwd.de 01326 ! 01327 ! History: 01328 ! Version Date Name 01329 ! ---------- ---------- ---- 01330 ! 1.00 2005/11/17 Dmitrii Mironov 01331 ! Initial release 01332 ! !VERSION! !DATE! <Your name> 01333 ! <Modification comments> 01334 ! 01335 ! Code Description: 01336 ! Language: Fortran 90. 01337 ! Software Standards: "European Standards for Writing and 01338 ! Documenting Exchangeable Fortran 90 Code". 01339 !============================================================================== 01340 ! 01341 ! Declarations: 01342 ! 01343 ! Modules used: 01344 01345 !_dm Parameters are USEd in module "sfcflx". 01346 !_nu USE modd_data_parameters , ONLY : & 01347 !_nu ireals, & ! KIND-type parameter for real variables 01348 !_nu iintegers ! KIND-type parameter for "normal" integer variables 01349 01350 !============================================================================== 01351 01352 IMPLICIT NONE 01353 01354 !============================================================================== 01355 ! 01356 ! Declarations 01357 01358 ! Input (function argument) 01359 REAL , INTENT(IN) :: 01360 T_dry , ! Dry air temperature [K] 01361 T_wetbulb , ! Wet bulb temperature [K] 01362 satwvpres_bulb , ! Satururation vapour pressure at wet-bulb temperature [N m^{-2}] 01363 P ! Atmospheric pressure [N m^{-2}] 01364 01365 ! Output (function result) 01366 REAL :: 01367 sfcflx_wvpreswetbulb ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 01368 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01369 01370 !============================================================================== 01371 ! Start calculations 01372 !------------------------------------------------------------------------------ 01373 01374 ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}] 01375 01376 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_WVPRESWETBULB',0,ZHOOK_HANDLE) 01377 sfcflx_wvpreswetbulb = satwvpres_bulb & 01378 - tpsf_c_a_p*P/tpsf_L_evap/tpsf_Rd_o_Rv*(T_dry-T_wetbulb) 01379 IF (LHOOK) CALL DR_HOOK('SFCFLX:SFCFLX_WVPRESWETBULB',1,ZHOOK_HANDLE) 01380 01381 01382 !------------------------------------------------------------------------------ 01383 ! End calculations 01384 !============================================================================== 01385 01386 END FUNCTION sfcflx_wvpreswetbulb 01387 01388 01389 END MODULE mode_sfcflx 01390
1.8.0
