File indexing completed on 2018-03-02 18:37:38 UTC
view on githubraw file Latest commit 672ec4d5 on 2016-12-01 17:05:48 UTC
672ec4d5cf Jean*0001 <HTML>
0002 <TITLE>module monin_obukhov</TITLE>
0003 <BODY BGCOLOR="#AABBCC" TEXT="#332211" >
0004
0005 <DIV ALIGN="CENTER"> <FONT SIZE=1>
0006 <A HREF="#INTERFACE">PUBLIC INTERFACE</A> /
0007 <A HREF="#ROUTINES">ROUTINES</A> /
0008 <A HREF="#NAMELIST">NAMELIST</A> /
0009 <A HREF="#CHANGES">CHANGES</A> /
0010 <A HREF="#ERRORS">ERRORS</A> /
0011 <A HREF="#REFERENCES">REFERENCES</A> /
0012 <A HREF="#NOTES">NOTES</A>
0013 </FONT>
0014 <BR><BR></DIV><HR>
0015
0016
0017 <H2>module monin_obukhov</H2>
0018 <A NAME="HEADER">
0019 <PRE>
0020 <B>Contact:</B> Isaac Held
0021 <B>Reviewers:</B>
0022
0023 <B><A HREF=".doc.log#monin_obukhov.f90">Tags/Status</A></B>
0024 </PRE>
0025 </A>
0026
0027 <A NAME="OVERVIEW">
0028 <HR>
0029 <H4>OVERVIEW</H4>
0030
0031 <PRE>
0032
0033 Routines for computing surface drag coefficients from
0034 data at the lowest model level using Monin-Obukhov scaling,
0035 for estimating the wind, temperature, and buoyancy profiles
0036 between the lowest model level and surface, and for computing
0037 diffusivities consistent these profiles
0038
0039 </PRE>
0040 </A>
0041
0042 <A NAME="DESCRIPTION">
0043
0044 <PRE>
0045
0046 Monin-Obukhov similarity (MOS) theory is the standard method
0047 for computing surface fluxes from the lowest level winds,
0048 temperatures, and tracer mixing ratios in GCMs. The lowest level
0049 is assumed to lie within the "surface layer" in which turbulent
0050 fluxes have negligible vertical variation, and in which MOS assumes
0051 that the wind and buoyancy profiles are a function only of the
0052 surface stress, the surface buoyancy flux, and the height z.
0053 A good reference is Garratt, J. R. "The Atmospheric Boundary Layer",
0054 Cambridge University Press, 1992. For a discussion of the detailed
0055 form of the similarity theory used in this module, see
0056 <A HREF="monin_obukhov.tech.ps">monin_obukhov.tech.ps</A>
0057
0058 The particular form of the monin-obukhov similarity theory utilized
0059 is defined by the similarity functions phi_m and phi_t chosen, where
0060 du/dzeta = phi_m * u_star/( vonkarm * zeta )
0061 db/dzeta( = phi_t * b_star/( vonkarm * zeta )
0062 (where b = buoyancy or specific humidity,
0063 vonkarm = Von Karman's constant
0064 u_star = friction velocity (stress = density*u_star**2)
0065 b_star = buoyancy scale (flux = density*u_star*b_star)
0066 zeta = z/L
0067 z = height of lowest model level
0068 L = monin_obukhov length )
0069 We use:
0070 on the unstable side
0071 phi_m = (1 - 16.0*zeta)**(-0.25)
0072 phi_t = (1 - 16.0*zeta)**(-0.5)
0073
0074 on the stable side -- phi_t = phi_m
0075 there are two options
0076
0077 option_1: phi = 1.0 + zeta*(5.0 + b_stab*zeta)/(1.0 + zeta)
0078 where b_stab = 1/rich_crit
0079 (rich_crit is a namelist parameter)
0080
0081 option_2: phi = 1 + 5.0*zeta (zeta < zeta_trans)
0082 = 1.0 + (5.0 - b_stab)*zeta_trans + b_stab*zeta
0083 (zeta_trans also a namelist parameter)
0084
0085 In order to change these similarity functions, one would need to
0086 1) modify the defintion of phi_m and phi_t in the internal
0087 subroutines mo_derivative_m and mo_derivative_t
0088 2) provide analytical versions of the integral similarity
0089 functions (PHI_m = int(zeta_0 to zeta) of phi/zeta)
0090 (see notes) in the subroutines mo_integral_m and mo_integral_tq
0091 3) modify stable_mix to be consistent with the new similarity functions
0092 (see notes) (mod_diff will take care of itself)
0093
0094
0095 </PRE>
0096 </A><!-- END DESCRIPTION -->
0097 <!--------------------------------------------------------------------->
0098 <A NAME="MODULES_USED">
0099 <HR>
0100 <H4>OTHER MODULES USED</H4>
0101 <!-- BEGIN MODULES_USED -->
0102 <PRE>
0103
0104 This modules uses
0105 constants_mod, only : grav, vonkarm
0106 utilities_mod, only : error_mesg, FATAL, file_exist, &
0107 check_nml_error, open_file, &
0108 get_my_pe, close_file
0109
0110
0111 </PRE>
0112 </A><!-- END MODULES_USED -->
0113 <!--------------------------------------------------------------------->
0114 <A NAME="INTERFACE">
0115 <HR>
0116 <H4>PUBLIC INTERFACE</H4>
0117 <!-- BEGIN INTERFACE -->
0118 <PRE>
0119
0120 use monin_obukhov_mod [,only: mo_drag , &
0121 mo_profile , &
0122 mo_diff , &
0123 stable_mix ]
0124
0125 mo_drag : computes the drag coefficients for momentum, heat, and moisture;
0126 also returns u_star (the friction velocity) and
0127 b_star (the buoyancy scale)
0128
0129 mo_profile : provides the profile of momentum, temperature, and
0130 specific humidity in the constant flux layer
0131 -- useful when one needs
0132 model output at a standard level below the lowest
0133 model level. For example, the term "surface wind"
0134 in meteorology often refers to
0135 the wind 10 meters above the surface.
0136
0137 mo_diff : computes the difusivity that, acting in isolation,
0138 would produce the correct profiles in ths surface layer
0139
0140 stable_mix : returns f(Ri) under stable conditions where Ri is the
0141 local Richardson's number, from which on can compute the
0142 diffusivity that produces the correct profile in the constant
0143 flux layer from D = l**2 |dv\dz| F(Ri) where l = vonkarm*z .
0144 Useful for matching a local Ri-dependent diffusivity in the free
0145 atmosphere under stable conditions o that implied by the
0146 surface layer profile
0147
0148
0149
0150 Notes:
0151
0152 all of the interfaces can be called with either 2d, 1d or 0d in/output.
0153 The base routines for mo_drag and mo_profile are 1d, while
0154 those for mo_diff and stable_mix are 2d, as these are the
0155 versions called within our GCMs. (mo_drag and mo_profile are
0156 called from the exchange_grid, which is 1d). Other versions call these
0157 base routines -- they are included for convenience -- the presumption is
0158 that efficiency is not of particular importance for these alternative
0159 versions)
0160
0161
0162
0163
0164 </PRE>
0165 </A>
0166
0167 <A NAME="ROUTINES">
0168 <HR>
0169 <H4>PUBLIC ROUTINES</H4>
0170
0171 <PRE>
0172
0173 ==========================================================================
0174
0175 call mo_drag (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, &
0176 u_star, b_star, [mask])
0177
0178 input:
0179
0180 pt, real, dimension(:)
0181 virtual potential temperature at lowest model level
0182 degrees Kelvin
0183
0184 pt0, real, dimension(:)
0185 virtual potential temperature at surface
0186 degrees Kelvin
0187
0188 z, real, dimension(:)
0189 height above surface of lowest model layer
0190 meters
0191
0192 z0, real, dimension(:)
0193 surface roughness for momentum
0194 meters
0195
0196 zt, real, dimension(:)
0197 surface roughness for temperature
0198 meters
0199
0200 zq, real, dimension(:)
0201 surface roughness for moisture
0202 meters
0203
0204 speed, real, dimension(:)
0205 wind speed at lowest model level with respect to surface
0206 (any "gustiness" factor should be included in speed)
0207 meters/sec
0208
0209 output:
0210
0211 drag_m, real, dimension(:)
0212 non-dimensional drag coefficient for momentum
0213
0214 drag_t, real, dimension(:)
0215 non-dimensional drag coefficient for temperature
0216
0217 drag_q, real, dimension(:)
0218 non-dimensional drag coefficient for specific humidity
0219
0220 u_star, real, dimension(:)
0221 friction velocity
0222 meters/sec
0223
0224 b_star, real, dimension(:)
0225 buoyancy scale
0226 (meters/sec)**2
0227
0228
0229 The magnitude of the wind stress is
0230 density*(ustar**2)
0231 The drag coefficient for momentum is
0232 (u_star/speed)**2
0233 The buoyancy flux is
0234 density*ustar*bstar
0235 The drag coefficient for heat etc is
0236 (u_star/speed)*(b_star/delta_b)
0237 where delta_b is the buoyancy difference between
0238 the surface and the lowest model level
0239 The buoyancy == grav*(p0 - p)/p0
0240
0241
0242 optional:
0243 mask : logical, dimension(:)
0244 computation performed only where mask = .true.
0245
0246
0247 Also:
0248
0249 can be called with all 1d arrays replaced by 2d arrays or by 0d scalars
0250 NOTE: in the 0d and 2d cases, there is no avail optional argument
0251 ==========================================================================
0252
0253 subroutine mo_profile(zref_m, zref_t, z, z0, zt, zq, &
0254 u_star, b_star, q_star, del_m, del_t, del_q, [mask])
0255
0256 input:
0257
0258 zref_m, real
0259 height above surface to which interpolation is requested
0260 for horizontal components of momentum
0261 meters
0262
0263 zref_t, real
0264 height above surface to which interpolation is requested
0265 for temperature and specific humidity
0266 meters
0267
0268 z, real, dimension(:)
0269 height of lowest model layer above surface
0270 meters
0271
0272 z0, real, dimension(:)
0273 surface roughness for momentum
0274 meters
0275
0276 zt, real, dimension(:)
0277 surface roughness for temperature
0278 meters
0279
0280 zq, real, dimension(:)
0281 surface roughness for moisture
0282 meters
0283
0284 u_star, real, dimension(:)
0285 friction velocity
0286 meters/sec
0287
0288 b_star, real, dimension(:)
0289 buoyancy scale
0290 (meters/sec)**2
0291
0292 q_star, real, dimension(:)
0293 moisture scale
0294 kg/kg
0295
0296 (Note: u_star and b_star are output from mo_drag,
0297 q_star = flux_q/(u_star * density)
0298
0299 optional input:
0300
0301 mask, logical, dimension(:)
0302 computation performed only where mask = .true.
0303
0304
0305 output:
0306
0307 del_m, real, dimension(:)
0308 dimensionless ratio, as defined below, for momentum
0309
0310 del_t, real, dimension(:)
0311 dimensionless ratio, as defined below, for temperature
0312
0313 del_q, real, dimension(:) o
0314 dimensionless ratio, as defined below, for specific humidity
0315
0316 Ratios are (f(zref) - f_surf)/(f(z) - f_surf)
0317
0318 Also:
0319
0320 can be called with all 1d arrays replaced by 2d arrays or by 0d scalars
0321 NOTE: in the 0d and 2d cases, there is no avail optional argument
0322
0323 in addition, this routine can be called with the interface
0324
0325 subroutine mo_profile(zref, z, z0, zt, zq, &
0326 u_star, b_star, q_star, del_m, del_h, del_q)
0327
0328 in which zref is a 1d array of the heights at which the profiles
0329 are requested, and all other input is 0d, 1d, or 2d as above,
0330 while the output del_m, del_h, del_q is dimensioned
0331 in the 0d case: dimension (size(zref))
0332 in the 1d case: dimension (size(input), size(zref))
0333 in the 2d case: dimension (size(input,1), size(input,2), size(zref))
0334 (note that the two inputs zref_m and zref_h and here replaced by
0335 the input vector zref)
0336
0337 ==========================================================================
0338
0339 subroutine mo_diff(z, u_star, b_star, k_m, k_h)
0340
0341 input:
0342
0343 u_star, real, dimension(:,:) -- 2d horizontal array
0344 surface friction velocity
0345 (meters/sec)
0346
0347 b_star, real, dimension(:,:) -- 2d horizontal array
0348 buoyancy scale
0349 (meters/sec**2)
0350
0351 (Note: u_star and b_star are output from mo_drag)
0352
0353 z, real, dimension(:,:,:)
0354 (first two dimensions conform to u_star, b_star)
0355 height above surface at which diffusivities
0356 are desired -- i.e., diffusivities returned at
0357 z(i,j,:) for each point (i,j)
0358 (meters)
0359
0360 output:
0361
0362 k_m : real, dimension(:,:,:) conforms to z
0363 kinematic diffusivity for momentum
0364 (meters**2/sec)
0365
0366 k_h : real, dimension(:,:,:)
0367 kinematic diffusivity for temperature
0368 (meters**2/sec)
0369
0370 Also 0d and 1d versions available
0371 0d: u_star, b_star = scalars;
0372 z, k_m, k_h = dimension(:)
0373 1d u_star, b_star = dimension(:),
0374 z, k_m, k_h = dimension(size(u_star), :)
0375
0376 Also, can be called so as to return diffusivities at one level
0377 (the height of this level can still vary from column to column)
0378 0d: u_star, b_star, z, k_m, k_h = scalars
0379 1d: u_star, b_star, z, k_m, k_h = dimension(:)
0380 2d u_star, b_star, z, k_m, k_h = dimension(:,:)
0381
0382 ==========================================================================
0383
0384 subroutine stable_mix(rich, mix)
0385
0386 input:
0387
0388 rich, real, dimension(:,:,:)
0389 local Richardson's number
0390 (none)
0391
0392 output:
0393
0394 mix, real, dimension(:,:,:)
0395 dimensional factor that produces the diffusivity
0396 consistent with the similarity profile
0397 under stable conditions when multiplied by
0398 (L**2) * |dv/dz| where v is the vector horizontal wind
0399 and L = vonkarm*z
0400
0401 Also 0d, 1d, and 2d versions vailable
0402
0403 </PRE>
0404 </A><!-- END ROUTINES -->
0405 <!--------------------------------------------------------------------->
0406 <A NAME="NAMELIST">
0407 <HR>
0408 <H4>NAMELIST</H4>
0409 <!-- BEGIN NAMELIST -->
0410 <PRE>
0411
0412 (default values provided)
0413
0414 logical :: neutral = .false.
0415
0416 If set to .true., all stability dependence is suppressed and
0417 neutral logarithmic profiles are used for all calculations;
0418 all other namelist parameters ignored
0419
0420 real :: drag_min = 1.e-05 (non-dimensional)
0421
0422 The drag coefficients (for both momentum and temperature)
0423 are not allowed to fall below drag_min
0424
0425
0426 integer :: stable_option = 1
0427
0428 must equal either 1 or 2
0429 two version for the similarity functions on the stable side
0430
0431 real :: rich_crit = 2.0 (non-dimensional)
0432
0433 The first step in applying the similarity theory is computing
0434 the bulk Richardson's number Ri between the lowest model level
0435 and the surface. If Ri is greater than rich_crit, the drag
0436 coefficients are set to drag_min.
0437
0438 This value also affects the drag coefficients under stable conditions
0439 for either choice of stable_option, as the drag is arranged to
0440 becomes very small as rich_crit is approached
0441
0442 real :: zeta_trans = 0.5 (non-dimensional)
0443
0444 A parameter in the stable-option = 2 piecewise linear similarity
0445 function for the stable side. Increasing zeta_trans reduces
0446 the drag and the implied diffusivities
0447
0448
0449 </PRE>
0450 </A>
0451
0452 <A NAME="CHANGES">
0453 <HR>
0454 <H4>CHANGE HISTORY</H4>
0455
0456 <PRE>
0457 <B><A HREF=".doc.log#monin_obukhov.f90">Revision history</A></B>
0458
0459 <b>changes</b> (10/4/1999)
0460
0461 MPP version created. Minor changes for open_file, error_mesg,
0462 and Fortran write statements. Answers should reproduce the
0463 previous version.
0464
0465 <b>more changes</b> (10/01)
0466
0467 -- time smoothing option removed
0468 -- second optional similarity function on stable side added
0469 -- stable_mix interface added
0470 -- code rearranged substantially for clarity, removing initial
0471 aggregation of points into unstable and stable sets
0472 -- answers reproduce for stable_option = 1 except in the very
0473 rare instance of points that are very close to neutral.
0474 These need to be treated separately to avoid a divide by sero.
0475 In this version, the drag coefficients at these points are simply
0476 set to the neutral value
0477
0478 </PRE>
0479 </A>
0480
0481 <A NAME="ERRORS">
0482 <HR>
0483 <H4>ERROR MESSAGES</H4>
0484
0485 <PRE>
0486
0487 FATAL ERRORS in MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD
0488
0489 -- rich_crit in monin_obukhov_mod must be > 0.25
0490
0491 -- drag_min in monin_obukhov_mod must be >= 0.0
0492
0493 -- the only allowable values of stable_option are 1 and 2
0494
0495 -- zeta_trans must be positive
0496
0497
0498 FATAL ERROR in solve_zeta in monin_obukhov_mod
0499
0500 -- surface drag iteration did not converge
0501
0502 </PRE>
0503 </A>
0504
0505 <A NAME="REFERENCES">
0506 <HR>
0507 <H4>REFERENCES</H4>
0508
0509 <PRE>
0510
0511 </PRE>
0512 </A>
0513
0514 <A NAME="BUGS">
0515 <HR>
0516 <H4>KNOWN BUGS</H4>
0517
0518 <PRE>
0519 None known
0520 </PRE>
0521 </A>
0522
0523 <A NAME="NOTES">
0524 <HR>
0525 <H4>NOTES</H4>
0526
0527 <PRE>
0528
0529 If iteration convergence is ever a problem, one can consider
0530 precomputation and table interpolation
0531
0532 It would be convenient if changing the derivative similarity functions
0533 automatically modified the integral functions and stable_mix
0534 (this would presumably require numerical integration and table look-ups
0535 as opposed to analytical expressions)
0536
0537
0538 </PRE>
0539 </A>
0540
0541 <A NAME="PLANS">
0542 <HR>
0543 <H4>FUTURE PLANS</H4>
0544
0545 <PRE>
0546
0547
0548 </PRE>
0549 </A>
0550
0551
0552 <HR>
0553 </BODY>
0554 </HTML>
