| blob | ac920798aca70820a87778bb363e0330d6d4e493 |
1 !WRF:MEDIATION_LAYER:PHYSICS
3 #if (NMM_CORE == 1)
4 MODULE diag_functions
5 CONTAINS
6 SUBROUTINE diag_functions_stub
7 END SUBROUTINE diag_functions_stub
8 END MODULE diag_functions
9 #else
11 MODULE diag_functions
13 CONTAINS
15 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
16 !~
17 !~ Name:
18 !~ calc_rh
19 !~
20 !~ Description:
21 !~ This function calculates relative humidity given pressure,
22 !~ temperature, and water vapor mixing ratio.
23 !~
24 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
25 FUNCTION calc_rh ( p, t, qv ) result ( rh )
27 IMPLICIT NONE
29 REAL, INTENT(IN) :: p, t, qv
30 REAL :: rh
32 ! Local
33 ! -----
34 REAL, PARAMETER :: pq0=379.90516
35 REAL, PARAMETER :: a2=17.2693882
36 REAL, PARAMETER :: a3=273.16
37 REAL, PARAMETER :: a4=35.86
38 REAL, PARAMETER :: rhmin=1.
39 REAL :: q, qs
40 INTEGER :: i,j,k
42 ! Following algorithms adapted from WRFPOST
43 ! May want to substitute with another later
44 ! -----------------------------------------
45 q=qv/(1.0+qv)
46 qs=pq0/p*exp(a2*(t-a3)/(t-a4))
47 rh=100.*q/qs
48 IF (rh .gt. 100.) THEN
49 rh=100.
50 ELSE IF (rh .lt. rhmin) THEN
51 rh=rhmin
52 ENDIF
54 END FUNCTION calc_rh
56 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
57 !~
58 !~ Name:
59 !~ uv_wind
60 !~
61 !~ Description:
62 !~ This function calculates the wind speed given U and V wind
63 !~ components.
64 !~
65 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
66 FUNCTION uv_wind ( u, v ) result ( wind_speed )
68 IMPLICIT NONE
70 REAL, INTENT(IN) :: u, v
71 REAL :: wind_speed
73 wind_speed = sqrt( u*u + v*v )
75 END FUNCTION uv_wind
79 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
80 !~
81 !~ Name:
82 !~ Theta
83 !~
84 !~ Description:
85 !~ This function calculates potential temperature as defined by
86 !~ Poisson's equation, given temperature and pressure ( hPa ).
87 !~
88 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
89 FUNCTION Theta ( t, p )
90 IMPLICIT NONE
92 !~ Variable declaration
93 ! --------------------
94 REAL, INTENT ( IN ) :: t
95 REAL, INTENT ( IN ) :: p
96 REAL :: theta
98 REAL :: Rd ! Dry gas constant
99 REAL :: Cp ! Specific heat of dry air at constant pressure
100 REAL :: p0 ! Standard pressure ( 1000 hPa )
102 Rd = 287.04
103 Cp = 1004.67
104 p0 = 1000.00
106 !~ Poisson's equation
107 ! ------------------
108 theta = t * ( (p0/p)**(Rd/Cp) )
110 END FUNCTION Theta
112 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
113 !~
114 !~ Name:
115 !~ Thetae
116 !~
117 !~ Description:
118 !~ This function returns equivalent potential temperature using the
119 !~ method described in Bolton 1980, Monthly Weather Review, equation 43.
120 !~
121 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
122 FUNCTION Thetae ( tK, p, rh, mixr )
123 IMPLICIT NONE
125 !~ Variable Declarations
126 ! ---------------------
127 REAL :: tK ! Temperature ( K )
128 REAL :: p ! Pressure ( hPa )
129 REAL :: rh ! Relative humidity
130 REAL :: mixr ! Mixing Ratio ( kg kg^-1)
131 REAL :: te ! Equivalent temperature ( K )
132 REAL :: thetae ! Equivalent potential temperature
134 REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
135 REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
136 REAL, PARAMETER :: lv = 2.54*(10**6) ! Latent heat of vaporization
137 ! (J kg^-1)
138 REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
139 ! at pressure (J/deg kg)
140 REAL :: tlc ! LCL temperature
142 !~ Calculate the temperature of the LCL
143 ! ------------------------------------
144 tlc = TLCL ( tK, rh )
146 !~ Calculate theta-e
147 ! -----------------
148 thetae = (tK * (p0/p)**( (R/Cp)*(1.- ( (.28E-3)*mixr*1000.) ) ) )* &
149 exp( (((3.376/tlc)-.00254))*&
150 (mixr*1000.*(1.+(.81E-3)*mixr*1000.)) )
152 END FUNCTION Thetae
156 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
157 !~
158 !~ Name:
159 !~ The2T.f90
160 !~
161 !~ Description:
162 !~ This function returns the temperature at any pressure level along a
163 !~ saturation adiabat by iteratively solving for it from the parcel
164 !~ thetae.
165 !~
166 !~ Dependencies:
167 !~ function thetae.f90
168 !~
169 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
170 FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
171 IMPLICIT NONE
173 !~ Variable Declaration
174 ! --------------------
175 REAL, INTENT ( IN ) :: thetaeK
176 REAL, INTENT ( IN ) :: pres
177 LOGICAL, INTENT ( INOUT ) :: flag
178 REAL :: tparcel
180 REAL :: thetaK
181 REAL :: tovtheta
182 REAL :: tcheck
183 REAL :: svpr, svpr2
184 REAL :: smixr, smixr2
185 REAL :: thetae_check, thetae_check2
186 REAL :: tguess_2, correction
188 LOGICAL :: found
189 INTEGER :: iter
191 REAL :: R ! Dry gas constant
192 REAL :: Cp ! Specific heat for dry air
193 REAL :: kappa ! Rd / Cp
194 REAL :: Lv ! Latent heat of vaporization at 0 deg. C
196 R = 287.04
197 Cp = 1004.67
198 Kappa = R/Cp
199 Lv = 2.500E+6
201 !~ Make initial guess for temperature of the parcel
202 ! ------------------------------------------------
203 tovtheta = (pres/100000.0)**(r/cp)
204 tparcel = thetaeK/exp(lv*.012/(cp*295.))*tovtheta
206 iter = 1
207 found = .false.
208 flag = .false.
210 DO
211 IF ( iter > 105 ) EXIT
213 tguess_2 = tparcel + REAL ( 1 )
215 svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
216 smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
217 svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
218 smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
220 ! ------------------------------------------------------------------ ~!
221 !~ When this function was orinially written, the final parcel ~!
222 !~ temperature check was based off of the parcel temperature and ~!
223 !~ not the theta-e it produced. As there are multiple temperature- ~!
224 !~ mixing ratio combinations that can produce a single theta-e value, ~!
225 !~ we change the check to be based off of the resultant theta-e ~!
226 !~ value. This seems to be the most accurate way of backing out ~!
227 !~ temperature from theta-e. ~!
228 !~ ~!
229 !~ Rentschler, April 2010 ~!
230 ! ------------------------------------------------------------------ !
232 !~ Old way...
233 !thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
234 !tcheck = thetaK * tovtheta
236 !~ New way
237 thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
238 thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
240 !~ Whew doggies - that there is some accuracy...
241 !IF ( ABS (tparcel-tcheck) < .05) THEN
242 IF ( ABS (thetaeK-thetae_check) < .001) THEN
243 found = .true.
244 flag = .true.
245 EXIT
246 END IF
248 !~ Old
249 !tparcel = tparcel + (tcheck - tparcel)*.3
251 !~ New
252 correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
253 tparcel = tparcel + correction
255 iter = iter + 1
256 END DO
258 !IF ( .not. found ) THEN
259 ! print*, "Warning! Thetae to temperature calculation did not converge!"
260 ! print*, "Thetae ", thetaeK, "Pressure ", pres
261 !END IF
263 END FUNCTION The2T
267 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
268 !~
269 !~ Name:
270 !~ VirtualTemperature
271 !~
272 !~ Description:
273 !~ This function returns virtual temperature given temperature ( K )
274 !~ and mixing ratio.
275 !~
276 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
277 FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
278 IMPLICIT NONE
280 !~ Variable declaration
281 real, intent ( in ) :: tK !~ Temperature
282 real, intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
283 real :: Tv !~ Virtual temperature
285 Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
287 END FUNCTION VirtualTemperature
292 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
293 !~
294 !~ Name:
295 !~ SaturationMixingRatio
296 !~
297 !~ Description:
298 !~ This function calculates saturation mixing ratio given the
299 !~ temperature ( K ) and the ambient pressure ( Pa ). Uses
300 !~ approximation of saturation vapor pressure.
301 !~
302 !~ References:
303 !~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
304 !~
305 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
306 FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
308 IMPLICIT NONE
310 REAL, INTENT ( IN ) :: tK
311 REAL, INTENT ( IN ) :: p
312 REAL :: ws
314 REAL :: es
316 es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
317 ws = ( 0.622*es ) / ( (p/100.0)-es )
319 END FUNCTION SaturationMixingRatio
323 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
324 !~
325 !~ Name:
326 !~ tlcl
327 !~
328 !~ Description:
329 !~ This function calculates the temperature of a parcel of air would have
330 !~ if lifed dry adiabatically to it's lifting condensation level (lcl).
331 !~
332 !~ References:
333 !~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
334 !~
335 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
336 FUNCTION TLCL ( tk, rh )
338 IMPLICIT NONE
340 REAL, INTENT ( IN ) :: tK !~ Temperature ( K )
341 REAL, INTENT ( IN ) :: rh !~ Relative Humidity ( % )
342 REAL :: tlcl
344 REAL :: denom, term1, term2
346 term1 = 1.0 / ( tK - 55.0 )
347 IF ( rh > REAL (0) ) THEN
348 term2 = ( LOG (rh/100.0) / 2840.0 )
349 ELSE
350 term2 = ( LOG (0.001/1.0) / 2840.0 )
351 END IF
352 denom = term1 - term2
353 tlcl = ( 1.0 / denom ) + REAL ( 55 )
355 END FUNCTION TLCL
359 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
360 !~
361 !~ Name:
362 !~ PWat
363 !~
364 !~ Description:
365 !~ This function calculates precipitable water by summing up the
366 !~ water vapor in a column from the first eta layer to model top
367 !~
368 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
369 FUNCTION Pwat ( nz, qv, qc, dz8w, rho )
371 IMPLICIT NONE
373 !~ Variable declaration
374 ! --------------------
375 INTEGER, INTENT ( IN ) :: nz !~ Number of vertical levels
376 REAL, INTENT ( IN ) :: qv ( nz ) !~ Specific humidity in layer (kg/kg)
377 REAL, INTENT ( IN ) :: qc ( nz ) !~ Cloud water in layer (kg/kg)
378 REAL, INTENT ( IN ) :: dz8w ( nz ) !~ Dist between vertical levels (m)
379 REAL, INTENT ( IN ) :: rho ( nz ) !~ Air density (kg/m^3)
380 REAL :: Pwat !~ Precipitable water (kg/m^2)
381 INTEGER :: k !~ Vertical index
383 !~ Precipitable water (kg/m^2)
384 ! ---------------------------
385 Pwat=0
386 DO k = 1, nz
387 !Based on AMS PWAT defination (https://glossary.ametsoc.org/wiki/Precipitable_water)
388 !PWAT is corrected as the column accumulated water vapor rather than water vapor + cloud water.
389 !Modified by Zhixiao
390 Pwat = Pwat + qv(k) * dz8w(k) * rho(k)
391 ENDDO
393 END FUNCTION Pwat
397 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
398 !~ ~!
399 !~ Name: ~!
400 !~ Buoyancy ~!
401 !~ ~!
402 !~ Description: ~!
403 !~ This function computes Convective Available Potential Energy (CAPE) ~!
404 !~ with inhibition as a result of water loading given the data required ~!
405 !~ to run up a sounding. ~!
406 !~ ~!
407 !~ Additionally, since we are running up a sounding anyways, this ~!
408 !~ function returns the height of the Level of Free Convection (LFC) and ~!
409 !~ the pressure at the LFC. That-a-ways, we don't have to run up a ~!
410 !~ sounding later, saving a relatively computationally expensive ~!
411 !~ routine. ~!
412 !~ ~!
413 !~ Usage: ~!
414 !~ ostat = Buoyancy ( tK, rh, p, hgt, sfc, CAPE, ZLFC, PLFC, parcel ) ~!
415 !~ ~!
416 !~ Where: ~!
417 !~ ~!
418 !~ IN ~!
419 !~ -- ~!
420 !~ tK = Temperature ( K ) ~!
421 !~ rh = Relative Humidity ( % ) ~!
422 !~ p = Pressure ( Pa ) ~!
423 !~ hgt = Geopotential heights ( m ) ~!
424 !~ sfc = integer rank within submitted arrays that represents the ~!
425 !~ surface ~!
426 !~ ~!
427 !~ OUT ~!
428 !~ --- ~!
429 !~ ostat INTEGER return status. Nonzero is bad. ~!
430 !~ CAPE ( J/kg ) Convective Available Potential Energy ~!
431 !~ ZLFC ( gpm ) Height at the LFC ~!
432 !~ PLFC ( Pa ) Pressure at the LFC ~!
433 !~ ~!
434 !~ tK, rh, p, and hgt are all REAL arrays, arranged from lower levels ~!
435 !~ to higher levels. ~!
436 !~ ~!
437 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
438 FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
439 parcel ) result (ostat)
441 IMPLICIT NONE
443 INTEGER, INTENT ( IN ) :: nz !~ Number of vertical levels
444 INTEGER, INTENT ( IN ) :: sfc !~ Surface level in the profile
445 REAL, INTENT ( IN ) :: tk ( nz ) !~ Temperature profile ( K )
446 REAL, INTENT ( IN ) :: rh ( nz ) !~ Relative Humidity profile ( % )
447 REAL, INTENT ( IN ) :: p ( nz ) !~ Pressure profile ( Pa )
448 REAL, INTENT ( IN ) :: hgt ( nz ) !~ Height profile ( gpm )
449 REAL, INTENT ( OUT ) :: cape !~ CAPE ( J kg^-1 )
450 REAL, INTENT ( OUT ) :: cin !~ CIN ( J kg^-1 )
451 REAL, INTENT ( OUT ) :: zlfc !~ LFC Height ( gpm )
452 REAL, INTENT ( OUT ) :: plfc !~ LFC Pressure ( Pa )
453 REAL, INTENT ( OUT ) :: lidx !~ Lifted index
454 INTEGER :: ostat !~ Function return status
455 !~ Nonzero is bad.
457 INTEGER, INTENT ( IN ) :: parcel !~ Most Unstable = 1 (default)
458 !~ Mean layer = 2
459 !~ Surface based = 3
461 !~ Derived profile variables
462 ! -------------------------
463 REAL :: ws ( nz ) !~ Saturation mixing ratio
464 REAL :: w ( nz ) !~ Mixing ratio
465 REAL :: etheta( nz )!~ Equivalent potential temperature. Modified by Zhixiao.
466 REAL :: dTvK ( nz ) !~ Parcel / ambient Tv difference
467 REAL :: buoy ( nz ) !~ Buoyancy
468 REAL :: tlclK !~ LCL temperature ( K )
469 REAL :: plcl !~ LCL pressure ( Pa )
470 REAL :: pel !~ Equilibrium pressure ( Pa ). Modified by Zhixiao.
471 REAL :: nbuoy !~ Negative buoyancy
472 REAL :: pbuoy !~ Positive buoyancy
474 !~ Source parcel information
475 ! -------------------------
476 REAL :: srctK !~ Source parcel temperature ( K )
477 REAL :: srcrh !~ Source parcel rh ( % )
478 REAL :: srcws !~ Source parcel sat. mixing ratio
479 REAL :: srcw !~ Source parcel mixing ratio
480 REAL :: srcp !~ Source parcel pressure ( Pa )
481 REAL :: srctheta !~ Source parcel theta ( K )
482 REAL :: srcthetaeK !~ Source parcel theta-e ( K )
483 INTEGER :: srclev !~ Level of the source parcel
484 REAL :: spdiff !~ Pressure difference
485 REAL :: srce !~ Equivalent potential temperature ( K ). Modified by Zhixiao.
487 !~ Parcel variables
488 ! ----------------
489 REAL :: ptK !~ Parcel temperature ( K )
490 REAL :: ptvK !~ Parcel virtual temperature ( K )
491 REAL :: tvK !~ Ambient virtual temperature ( K )
492 REAL :: pw !~ Parcel mixing ratio
494 !~ Other utility variables
495 ! -----------------------
496 INTEGER :: i, j, k !~ Dummy iterator
497 INTEGER :: lfclev !~ Level of LFC
498 INTEGER :: ellev !~ Level of EL. Modified by Zhixiao.
499 INTEGER :: prcl !~ Internal parcel type indicator
500 INTEGER :: mlev !~ Level for ML calculation
501 INTEGER :: lyrcnt !~ Number of layers in mean layer
502 LOGICAL :: flag !~ Dummy flag
503 LOGICAL :: wflag !~ Saturation flag
504 REAL :: freeze !~ Water loading multiplier
505 REAL :: pdiff !~ Pressure difference between levs
506 REAL :: pm, pu, pd !~ Middle, upper, lower pressures
507 REAL :: lidxu !~ Lifted index at upper level
508 REAL :: lidxd !~ Lifted index at lower level
510 !~ Thermo / dynamical constants
511 ! ----------------------------
512 REAL :: Rd !~ Dry gas constant
513 PARAMETER ( Rd = 287.058 ) !~ J deg^-1 kg^-1
514 REAL :: Cp !~ Specific heat constant pressure
515 PARAMETER ( Cp = 1004.67 ) !~ J deg^-1 kg^-1
516 REAL :: g !~ Acceleration due to gravity
517 PARAMETER ( g = 9.80665 ) !~ m s^-2
518 REAL :: RUNDEF
519 PARAMETER ( RUNDEF = -9.999E30 )
520 !~ Initialize variables
521 ! --------------------
522 ostat = 0
523 CAPE = REAL ( 0 )
524 CIN = RUNDEF !Change CIN filling values from 0 to default filling. CIN should not initially be filled by 0, because 0 means no inhibition energy. Modified by Zhixiao
525 ZLFC = RUNDEF
526 PLFC = RUNDEF
528 !~ Look for submitted parcel definition
529 !~ 1 = Most unstable
530 !~ 2 = Mean layer
531 !~ 3 = Surface based
532 ! -------------------------------------
533 IF ( parcel > 3 .or. parcel < 1 ) THEN
534 prcl = 1
535 ELSE
536 prcl = parcel
537 END IF
539 !~ Initalize our parcel to be (sort of) surface based. Because of
540 !~ issues we've been observing in the WRF model, specifically with
541 !~ excessive surface moisture values at the surface, using a true
542 !~ surface based parcel is resulting a more unstable environment
543 !~ than is actually occuring. To address this, our surface parcel
544 !~ is now going to be defined as the parcel between 25-50 hPa
545 !~ above the surface. UPDATE - now that this routine is in WRF,
546 !~ going to trust surface info. GAC 20140415
547 ! ----------------------------------------------------------------
549 !~ Compute mixing ratio values for the layer
550 ! -----------------------------------------
551 DO k = sfc, nz
552 ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
553 w ( k ) = ( rh(k)/100.0 ) * ws ( k )
554 !Removed by Zhixiao. Water vapor mixing ratio (w) is not conserved during parcel lifting processes. We should not use w to define MU layer.
555 !thetav(k) = Theta ( VirtualTemperature ( tK (k), w (k) ), p(k)/100.0 )
556 !Added by Zhixiao. Critical modification: We use the model level with maximum equivalent potential temperature (etheta) below 500mb to define the MU layer
557 !Because equivalent potential temperature is conserved in dry and moist adiabatic processes.
558 etheta(k) = Thetae( tK(k), p(k)/100.0, rh(k), w(k) )
559 END DO
561 srclev = sfc
562 srctK = tK ( sfc )
563 srcrh = rh ( sfc )
564 srcp = p ( sfc )
565 srcws = ws ( sfc )
566 srcw = w ( sfc )
567 srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
568 srce = etheta (sfc) !Modified by Zhixiao
570 !~ Compute the profile mixing ratio. If the parcel is the MU parcel,
571 !~ define our parcel to be the most unstable parcel within the lowest
572 !~ 180 mb.
573 ! -------------------------------------------------------------------
574 mlev = sfc + 1
575 !Change initial searching level from the second to first model level. Because we did not compute pdiff, and p(k-1) properties is unnecessary.
576 !Modified by Zhixiao.
577 DO k = sfc, nz
579 !~ Identify the last layer within 100 hPa of the surface
580 ! -----------------------------------------------------
581 pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
582 IF ( pdiff <= REAL (100) ) mlev = k
584 !~ If we've made it past the lowest 500 hPa, exit the loop. MU layer is assumed below 500 hPa. Modified by Zhixiao.
585 ! -------------------------------------------------------
586 IF ( p(k) <= REAL (50000) ) EXIT
588 IF ( prcl == 1 ) THEN
589 ! Removed by Zhixiao, w can not used for defining MU layer
590 !IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
591 ! Modified by Zhixiao, MU layer is featured by the max etheta
592 IF (etheta(k) > srce) THEN !Modified by Zhixiao.
593 srctheta = Theta ( tK(k), p(k)/100.0 )
594 srcw = w ( k )
595 srclev = k
596 srctK = tK ( k )
597 srcrh = rh ( k )
598 srcp = p ( k )
599 srce = etheta(k) !Modified by Zhixiao
600 END IF
601 END IF
603 END DO
605 !~ If we want the mean layer parcel, compute the mean values in the
606 !~ lowest 100 hPa.
607 ! ----------------------------------------------------------------
608 lyrcnt = mlev - sfc + 1
609 IF ( prcl == 2 ) THEN
611 srclev = sfc
612 srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
613 srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
614 srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
615 srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
616 srctheta = Theta ( srctK, srcp/100. )
618 END IF
620 srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
621 !~ Calculate temperature and pressure of the LCL
622 ! ---------------------------------------------
623 tlclK = TLCL ( tK(srclev), rh(srclev) )
624 plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
626 !~ Now lift the parcel
627 ! -------------------
629 buoy = REAL ( 0 )
630 pw = srcw
631 wflag = .false.
632 DO k = srclev, nz
633 IF ( p (k) <= plcl ) THEN
635 !~ The first level after we pass the LCL, we're still going to
636 !~ lift the parcel dry adiabatically, as we haven't added the
637 !~ the required code to switch between the dry adiabatic and moist
638 !~ adiabatic cooling. Since the dry version results in a greater
639 !~ temperature loss, doing that for the first step so we don't over
640 !~ guesstimate the instability.
641 ! ----------------------------------------------------------------
643 IF ( wflag ) THEN
644 flag = .false.
646 !~ Above the LCL, our parcel is now undergoing moist adiabatic
647 !~ cooling. Because of the latent heating being undergone as
648 !~ the parcel rises above the LFC, must iterative solve for the
649 !~ parcel temperature using equivalant potential temperature,
650 !~ which is conserved during both dry adiabatic and
651 !~ pseudoadiabatic displacements.
652 ! --------------------------------------------------------------
653 ptK = The2T ( srcthetaeK, p(k), flag )
655 !~ Calculate the parcel mixing ratio, which is now changing
656 !~ as we condense moisture out of the parcel, and is equivalent
657 !~ to the saturation mixing ratio, since we are, in theory, at
658 !~ saturation.
659 ! ------------------------------------------------------------
660 pw = SaturationMixingRatio ( ptK, p(k) )
662 !~ Now we can calculate the virtual temperature of the parcel
663 !~ and the surrounding environment to assess the buoyancy.
664 ! ----------------------------------------------------------
665 ptvK = VirtualTemperature ( ptK, pw )
666 tvK = VirtualTemperature ( tK (k), w (k) )
668 !~ Modification to account for water loading
669 ! -----------------------------------------
670 freeze = 0.033 * ( 263.15 - pTvK )
671 IF ( freeze > 1.0 ) freeze = 1.0
672 IF ( freeze < 0.0 ) freeze = 0.0
674 !~ Approximate how much of the water vapor has condensed out
675 !~ of the parcel at this level
676 ! ---------------------------------------------------------
677 freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
679 pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
680 dTvK ( k ) = ptvK - tvK
681 buoy ( k ) = g * ( dTvK ( k ) / tvK )
683 ELSE
685 !~ Since the theta remains constant whilst undergoing dry
686 !~ adiabatic processes, can back out the parcel temperature
687 !~ from potential temperature below the LCL
688 ! --------------------------------------------------------
689 ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
691 !~ Grab the parcel virtual temperture, can use the source
692 !~ mixing ratio since we are undergoing dry adiabatic cooling
693 ! ----------------------------------------------------------
694 ptvK = VirtualTemperature ( ptK, srcw )
696 !~ Virtual temperature of the environment
697 ! --------------------------------------
698 tvK = VirtualTemperature ( tK (k), w (k) )
700 !~ Buoyancy at this level
701 ! ----------------------
702 dTvK ( k ) = ptvK - tvK
703 buoy ( k ) = g * ( dtvK ( k ) / tvK )
705 wflag = .true.
707 END IF
709 ELSE
711 !~ Since the theta remains constant whilst undergoing dry
712 !~ adiabatic processes, can back out the parcel temperature
713 !~ from potential temperature below the LCL
714 ! --------------------------------------------------------
715 ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
717 !~ Grab the parcel virtual temperture, can use the source
718 !~ mixing ratio since we are undergoing dry adiabatic cooling
719 ! ----------------------------------------------------------
720 ptvK = VirtualTemperature ( ptK, srcw )
722 !~ Virtual temperature of the environment
723 ! --------------------------------------
724 tvK = VirtualTemperature ( tK (k), w (k) )
726 !~ Buoyancy at this level
727 ! ---------------------
728 dTvK ( k ) = ptvK - tvK
729 buoy ( k ) = g * ( dtvK ( k ) / tvK )
731 END IF
733 !~ Chirp
734 ! -----
735 ! WRITE ( *,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw*1000.,ptvK,tvK,buoy(k)
737 END DO
739 !~ Add up the buoyancies, find the LFC
740 ! -----------------------------------
741 flag = .false.
742 lfclev = -1
743 ellev = -1 !Modified by Zhixiao
744 DO k = sfc, nz !Modified by Zhixiao
745 !~ LFC is defiend as the highest level when negative buyancy turns postive.
746 ! -----------------------------------
747 IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) <= plcl ) THEN !Modified by Zhixiao
748 flag = .true.
749 lfclev = k
750 END IF
751 !~ Take the Highest EL as final result. Modified by Zhixiao
752 ! ----------------------------------------------------------------
753 IF (k >= 2) THEN !Modified by Zhixiao
754 IF ( flag .and. buoy (k) < REAL (0) .and. buoy (k-1) >= REAL (0)) THEN !Modified by Zhixiao
755 ellev = k !Modified by Zhixiao
756 END IF
757 END IF
758 ! When buoy turns negative again, reset LFC flag and keep the highest LFC as the effective output
759 IF (buoy (k) < REAL (0) .and. flag) THEN
760 flag = .false.
761 END IF
762 END DO
763 IF (ellev >= 0) THEN !Modified by Zhixiao
764 pel = p(ellev) !Modified by Zhixiao
765 CIN=REAL ( 0 ) !Modified by Zhixiao
766 DO k = sfc+1, nz
767 ! Make CAPE and CIN consistent with AMS definition
768 ! https://glossary.ametsoc.org/wiki/Convective_available_potential_energy
769 ! https://glossary.ametsoc.org/wiki/Convective_inhibition
770 IF ( p (k) <= plcl .and. p (k) > plfc) THEN !Modified by Zhixiao
771 ! CIN is the vertically integrated negative buoyant energy between LCL and LFC
772 CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
773 END IF
774 IF ( p (k) <= plfc .and. p (k) > pel) THEN !Modified by Zhixiao
775 ! CAPE is the vertically integrated positive buoyant energy between LFC and EL
776 CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
777 END IF !Modified by Zhixiao
778 END DO !Modified by Zhixiao
779 END IF !Modified by Zhixiao
780 !~ Calculate lifted index by interpolating difference between
781 !~ parcel and ambient Tv to 500mb.
782 ! ----------------------------------------------------------
783 DO k = sfc + 1, nz
785 pm = 50000.
786 pu = p ( k )
787 pd = p ( k - 1 )
789 !~ If we're already above 500mb just set lifted index to 0.
790 !~ --------------------------------------------------------
791 IF ( pd .le. pm ) THEN
792 lidx = 0.
793 EXIT
795 ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
797 !~ Found trapping pressure: up, middle, down.
798 !~ We are doing first order interpolation.
799 ! ------------------------------------------
800 lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
801 lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
802 lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
803 EXIT
805 ENDIF
807 END DO
808 !~ Assuming the the LFC is at a pressure level for now
809 ! ---------------------------------------------------
810 IF ( lfclev > 0 ) THEN
811 PLFC = p ( lfclev )
812 ZLFC = hgt ( lfclev )
813 END IF
815 IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN
816 PLFC = REAL ( -1 )
817 ZLFC = REAL ( -1 )
818 END IF
820 IF ( CAPE /= CAPE ) cape = REAL ( 0 )
822 IF ( CIN /= CIN ) cin = RUNDEF
824 !~ Chirp
825 ! -----
826 ! WRITE ( *,* ) ' CAPE: ', cape, ' CIN: ', cin
827 ! WRITE ( *,* ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
828 ! WRITE ( *,* ) ''
829 ! WRITE ( *,* ) ' Exiting buoyancy.'
830 ! WRITE ( *,* ) ' ==================================== '
831 ! WRITE ( *,* ) ''
833 END FUNCTION Buoyancy
837 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
838 ! . . .
839 ! SUBPROGRAM: NGMSLP NMC SEA LEVEL PRESSURE REDUCTION
840 ! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-02
841 !
842 ! ABSTRACT:
843 !
844 ! THIS ROUTINE COMPUTES SEA LEVEL PRESSURE USING THE
845 ! HYDROSTATIC EQUATION WITH THE SHUELL CORRECTION. THE
846 ! FOLLOWING IS BASED ON DOCUMENTATION IN SUBROUTINE
847 ! OUTHYDRO OF THE NGM:
848 !
849 ! THE FUNDAMENTAL HYDROSTATIC EQUATION IS
850 ! D(HEIGHT)
851 ! --------- = TAU = VIRTUAL TEMPERATURE * (RGAS/GRAVITY)
852 ! D (Z)
853 ! WHERE
854 ! Z = MINUS LOG OF PRESSURE (-LN(P)).
855 !
856 ! SEA-LEVEL PRESSURE IS COMPUTED FROM THE FORMULA
857 ! PRESS(MSL) = PRESS(GROUND) * EXP( F)
858 ! WHERE
859 ! F = HEIGHT OF GROUND / MEAN TAU
860 ! MEAN TAU = ( TAU(GRND) + TAU(SL) ) / 2
861 !
862 ! IN THE NGM TAU(GRND) AND TAU(SL) ARE FIRST SET USING A
863 ! 6.5DEG/KM LAPSE RATE FROM LOWEST MDL LEVEL. THIS IS MODIFIED
864 ! BY A CORRECTION BASED ON THE CRITICAL TAU OF THE SHUELL
865 ! CORRECTION:
866 ! TAUCR=(RGASD/GRAVITY) * 290.66
867 !
868 ! 1) WHERE ONLY TAU(SL) EXCEEDS TAUCR, CHANGE TAU(SL) TO TAUCR.
869 !
870 ! 2) WHERE BOTH TAU(SL) AND TAU(GRND) EXCEED TAUCR,
871 ! CHANGE TAU(SL) TO TAUCR-CONST*(TAU(GRND)-TAUCR )**2
872 ! WHERE CONST = .005 (GRAVITY/RGASD)
873 !
874 ! THE AVERAGE OF TAU(SL) AND TAU(GRND) IS THEN USED TOGETHER
875 ! WITH THE GROUND HEIGHT AND PRESSURE TO DERIVE THE PRESSURE
876 ! AT SEA LEVEL.
877 !
878 ! HEIGHT OF THE 1000MB SURFACE IS COMPUTED FROM THE MSL PRESSURE
879 ! FIELD USING THE FORMULA:
880 !
881 ! P(MSL) - P(1000MB) = MEAN DENSITY * GRAVITY * HGT(1000MBS)
882 !
883 ! WHERE P(MSL) IS THE SEA LEVEL PRESSURE FIELD WE HAVE JUST
884 ! COMPUTED.
885 !
886 !
887 ! MEB 6/13/02: THIS CODE HAS BEEN SIMPLIFIED CONSIDERABLY FROM
888 ! THE ONE USED IN ETAPOST. HORIZONTAL SMOOTHING HAS BEEN
889 ! REMOVED AND THE FIRST MODEL LEVEL IS USED RATHER
890 ! THAN THE MEAN OF THE VIRTUAL TEMPERATURES IN
891 ! THE LOWEST 30MB ABOVE GROUND TO COMPUTE TAU(GRND).
892 !
893 ! .
894 !
895 ! PROGRAM HISTORY LOG:
896 ! 93-02-02 RUSS TREADON
897 ! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D
898 ! 00-01-04 JIM TUCCILLO - MPI VERSION
899 ! 01-10-25 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
900 ! 01-11-02 H CHUANG - MODIFIED LINE 234 FOR COMPUTATION OF
901 ! SIGMA/HYBRID SLP
902 ! 01-12-18 H CHUANG - INCLUDED SMOOTHING ALONG BOUNDARIES TO BE
903 ! CONSISTENT WITH MESINGER SLP
904 ! 02-06-13 MIKE BALDWIN - WRF VERSION
905 ! 06-12-18 H CHUANG - BUG FIX TO CORRECT TAU AT SFC
906 ! 14-04-17 G CREIGHTON - MODIFIED TO INSERT INTO AFWA DIAGNOSTICS IN WRF
907 !
908 !$$$
910 FUNCTION MSLP ( zsfc, psfc, zlev1, qlev1, tlev1 )
912 implicit none
915 ! DECLARE VARIABLES
917 REAL, INTENT ( IN ) :: zsfc !~ Surface height ( m )
918 REAL, INTENT ( IN ) :: psfc !~ Surface height ( m )
919 REAL, INTENT ( IN ) :: zlev1 !~ Level 1 height ( m )
920 REAL, INTENT ( IN ) :: qlev1 !~ Level 1 mixing ratio ( kg/kg )
921 REAL, INTENT ( IN ) :: tlev1 !~ Level 1 temperature ( K )
922 real,PARAMETER :: G=9.81
923 real,PARAMETER :: GI=1./G
924 real,PARAMETER :: RD=287.0
925 real,PARAMETER :: ZSL=0.0
926 real,PARAMETER :: TAUCR=RD*GI*290.66,CONST=0.005*G/RD
927 real,PARAMETER :: GORD=G/RD,DP=60.E2
928 real,PARAMETER :: GAMMA=6.5E-3
930 real MSLP,TVRT,TVRSFC,TAUSFC,TVRSL,TAUSL,TAUAVG
931 !
932 !**********************************************************************
933 ! START NGMSLP HERE.
934 !
935 MSLP = PSFC
936 !
937 ! COMPUTE LAYER TAU (VIRTUAL TEMP*RD/G).
938 TVRT = TLEV1*(1.0+0.608*QLEV1)
939 !TAU = TVRT*RD*GI
940 !
941 ! COMPUTE TAU AT THE GROUND (Z=ZSFC) AND SEA LEVEL (Z=0)
942 ! ASSUMING A CONSTANT LAPSE RATE OF GAMMA=6.5DEG/KM.
943 TVRSFC = TVRT + (ZLEV1 - ZSFC)*GAMMA
944 TAUSFC = TVRSFC*RD*GI
945 TVRSL = TVRT + (ZLEV1 - ZSL)*GAMMA
946 TAUSL = TVRSL*RD*GI
947 !
948 ! IF NEED BE APPLY SHEULL CORRECTION.
949 IF ((TAUSL.GT.TAUCR).AND.(TAUSFC.LE.TAUCR)) THEN
950 TAUSL=TAUCR
951 ELSEIF ((TAUSL.GT.TAUCR).AND.(TAUSFC.GT.TAUCR)) THEN
952 TAUSL = TAUCR-CONST*(TAUSFC-TAUCR)**2
953 ENDIF
954 !
955 ! COMPUTE MEAN TAU.
956 TAUAVG = 0.5*(TAUSL+TAUSFC)
957 !
958 ! COMPUTE SEA LEVEL PRESSURE.
959 MSLP = PSFC*EXP(ZSFC/TAUAVG)
961 END FUNCTION MSLP
965 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
966 !~ ~!
967 !~ Name: ~!
968 !~ calc_fits ~!
969 !~ ~!
970 !~ Description: ~!
971 !~ This function computes Fighter Index Thermal Stress values given ~!
972 !~ dry bulb temperature, relative humidity, and pressure. ~!
973 !~ ~!
974 !~ Usage: ~!
975 !~ fitsval = calc_fits ( p, tK, rh ) ~!
976 !~ ~!
977 !~ Where: ~!
978 !~ p = Pressure ( Pa ) ~!
979 !~ tK = Temperature ( K ) ~!
980 !~ rh = Relative Humidity ( % ) ~!
981 !~ ~!
982 !~ Reference: ~!
983 !~ Stribley, R.F., S. Nunneley, 1978: Fighter Index of Thermal Stress: ~!
984 !~ Development of interim guidance for hot-weather USAF operations. ~!
985 !~ SAM-TR-78-6. Eqn. 9 ~!
986 !~ ~!
987 !~ Formula: ~!
988 !~ FITS = 0.8281*Twb + 0.3549*Tdb + 5.08 (degrees Celsius) ~!
989 !~ ~!
990 !~ Where: ~!
991 !~ Twb = Wet Bulb Temperature ~!
992 !~ Tdb = Dry Bulb Temperature ~!
993 !~ ~!
994 !~ Written: ~!
995 !~ Scott Rentschler, Software Engineering Services ~!
996 !~ Fine Scale Models Team ~!
997 !~ Air Force Weather Agency, 16WS/WXN ~!
998 !~ DSN: 271-3331 Comm: (402) 294-3331 ~!
999 !~ scott.rentschler@offutt.af.mil ~!
1000 !~ ~!
1001 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1002 FUNCTION calc_fits ( p, tK, rh ) RESULT ( fits )
1004 implicit none
1006 !~ Variable declaration
1007 ! --------------------
1008 real, intent ( in ) :: p !~ Pressure ( Pa )
1009 real, intent ( in ) :: tK !~ Temperature ( K )
1010 real, intent ( in ) :: rh !~ Rel Humidity ( % )
1011 real :: fits !~ FITS index value
1013 !~ Utility variables
1014 ! --------------------------
1015 real :: twb !~ Wet bulb temperature ( C )
1016 real :: wbt
1018 ! ---------------------------------------------------------------------- !
1019 ! ---------------------------------------------------------------------- !
1021 !~ Initialize variables
1022 ! --------------------
1023 fits = REAL ( 0 )
1025 !~ Get the wet bulb temperature in degrees Celsius
1026 ! -----------------------------------------------
1027 twb = WetBulbTemp ( p, tK, rh ) - 273.15
1029 !~ Compute the FITS index
1030 ! ----------------------
1031 fits = 0.8281*twb + 0.3549*( tK - 273.15 ) + 5.08
1033 !~ Convert the index to Kelvin
1034 ! ---------------------------
1035 fits = fits + 273.15
1037 END FUNCTION calc_fits
1041 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1042 !~
1043 !~ Name:
1044 !~ calc_wc
1045 !~
1046 !~ Description:
1047 !~ This function calculates wind chill given temperature ( K ) and
1048 !~ wind speed ( m/s )
1049 !~
1050 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1051 FUNCTION calc_wc ( tK, wspd ) RESULT ( wc )
1053 implicit none
1055 !~ Variable Declarations
1056 ! ---------------------
1057 real, intent ( in ) :: tK
1058 real, intent ( in ) :: wspd
1060 real :: tF, wc, wspd_mph
1062 wspd_mph = wspd * 2.23693629 ! convert to mph
1063 tF = (( tK - 273.15 ) * ( REAL (9) / REAL (5) ) ) + REAL ( 32 )
1065 wc = 35.74 &
1066 + ( 0.6215 * tF ) &
1067 - ( 35.75 * ( wspd_mph**0.16 ) ) &
1068 + ( 0.4275 * tF * ( wspd_mph**0.16 ) )
1070 wc = (( wc - REAL (32) ) * ( REAL (5) / REAL (9) ) ) + 273.15
1072 END FUNCTION calc_wc
1076 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1077 !~
1078 !~ Name:
1079 !~ calc_hi
1080 !~
1081 !~ Description:
1082 !~ This subroutine calculates the heat index. Requires temperature ( K )
1083 !~ and relative humidity ( % ).
1084 !~
1085 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1086 FUNCTION calc_hi ( Tk, RH ) result ( HI )
1088 implicit none
1090 !~ Variable declarations
1091 ! ---------------------
1092 real, intent ( in ) :: Tk
1093 real, intent ( in ) :: RH
1095 real :: tF, tF2, rh2, HI
1097 !~ If temperature > 70F then calculate heat index, else set it equal
1098 !~ to dry temperature
1099 ! -----------------------------------------------------------------
1100 IF ( Tk > 294.26111 ) THEN
1102 tF = ( (Tk - 273.15) * (REAL (9)/REAL (5)) ) + REAL ( 32 )
1103 tF2 = tF ** 2
1104 rh2 = RH ** 2
1106 HI = -42.379 &
1107 + ( 2.04901523 * tF ) &
1108 + ( 10.14333127 * RH ) &
1109 - ( 0.22475541 * tF * RH ) &
1110 - ( 6.83783E-03 * tF2 ) &
1111 - ( 5.481717E-02 * rh2 ) &
1112 + ( 1.22874E-03 * tF2 * RH ) &
1113 + ( 8.5282E-04 * tF * rh2 ) &
1114 - ( 1.99E-06 * tF2 * rh2 )
1116 HI = ((HI - REAL (32)) * (REAL (5)/REAL (9))) + 273.15
1117 ELSE
1118 HI = Tk
1119 END IF
1121 END FUNCTION calc_hi
1123 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1124 !~ ~!
1125 !~ Name: ~!
1126 !~ WetBulbTemp ~!
1127 !~ ~!
1128 !~ Description: ~!
1129 !~ This function approximates the Wet Bulb Temperature (K) provided ~!
1130 !~ dry bulb temperature (K), relative humidity (%), and pressure (Pa). ~!
1131 !~ ~!
1132 !~ Usage: ~!
1133 !~ wbt = WetBulbTemperature ( p, tK, rh ) ~!
1134 !~ ~!
1135 !~ Where: ~!
1136 !~ p = Pressure ( Pa ) ~!
1137 !~ tK = Temperature ( K ) ~!
1138 !~ rh = Relative Humidity ( % ) ~!
1139 !~ ~!
1140 !~ Reference: ~!
1141 !~ American Society of Civil Engineers ~!
1142 !~ Evapotraspiration and Irrigation Water Requirements ~!
1143 !~ Jensen et al (1990) ASCE Manual No. 70, pp 176-177 ~!
1144 !~ ~!
1145 !~ Written: ~!
1146 !~ Scott Rentschler, Software Engineering Services ~!
1147 !~ Fine Scale Models Team ~!
1148 !~ Air Force Weather Agency ~!
1149 !~ DSM: 271-3331 Comm: (402) 294-3331 ~!
1150 !~ scott.rentschler@offutt.af.mil ~!
1151 !~ ~!
1152 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1153 FUNCTION WetBulbTemp ( p, tK, rh) result( wbt )
1155 implicit none
1157 !~ Variable delclaration
1158 ! ---------------------
1159 real, intent ( in ) :: p !~ Pressure ( Pa )
1160 real, intent ( in ) :: tK !~ Temperature ( K )
1161 real, intent ( in ) :: rh !~ Relative Humidity ( % )
1162 real :: wbt !~ Wet Bulb Temperature ( K )
1164 !~ Utility variables
1165 ! -----------------
1166 real :: tdK !~ Dewpoint temperature ( K )
1167 real :: tC !~ Temperature ( C )
1168 real :: tdC !~ Dewpoint temperature ( K )
1169 real :: svapr !~ Saturation vapor pressure ( Pa )
1170 real :: vapr !~ Ambient vapor pressure ( Pa )
1171 real :: gamma !~ Dummy term
1172 real :: delta !~ Dummy term
1174 ! ---------------------------------------------------------------------- !
1175 ! ---------------------------------------------------------------------- !
1176 !~ Initialize variables
1177 ! --------------------
1178 wbt = REAL ( 0 )
1179 tC = tK - 273.15
1181 !~ Compute saturation vapor pressure ( Pa )
1182 ! ----------------------------------------
1183 svapr = calc_es ( tK ) * REAL ( 100 )
1185 !~ Compute vapor pressure
1186 ! ----------------------
1187 vapr = svapr * ( rh / REAL (100) )
1189 !~ Grab the dewpoint
1190 ! -----------------
1191 tdC = calc_Dewpoint ( tC, rh )
1192 tdK = tdC + 273.15
1194 !~ Compute dummy terms
1195 ! -------------------
1196 gamma = 0.00066 * ( p / REAL (1000) )
1197 delta = REAL ( 4098 ) * ( vapr / REAL(1000) ) / ( (tC+237.3)**2 )
1199 !~ Compute the wet bulb temperature
1200 ! --------------------------------
1201 wbt = ( ((gamma * tC) + (delta * tdC)) / (gamma + delta) ) + 273.15
1203 END FUNCTION WetBulbTemp
1206 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1207 !~
1208 !~ Name:
1209 !~ calc_Dewpoint
1210 !~
1211 !~ Description:
1212 !~ This function approximates dewpoint given temperature and rh.
1213 !~
1214 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1215 FUNCTION calc_Dewpoint ( tC, rh) result( Dewpoint )
1217 implicit none
1219 !~ Variable Declaration
1220 ! --------------------
1221 real, intent ( in ) :: tC
1222 real, intent ( in ) :: rh
1223 real :: Dewpoint
1225 real :: term, es, e1, e, logs, expon
1227 expon = ( 7.5*tC ) / ( 237.7+tC )
1228 es = 6.112 * ( 10**expon ) ! Saturated vapor pressure
1229 e = es * ( rh/100.0 ) ! Vapor pressure
1230 logs = LOG10 ( e/6.112 )
1231 Dewpoint = ( 237.7*logs ) / ( 7.5-logs )
1233 END FUNCTION calc_Dewpoint
1236 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1237 !~
1238 !~ Name:
1239 !~ calc_es
1240 !~
1241 !~ Description:
1242 !~ This function returns the saturation vapor pressure over water ( hPa )
1243 !~ given temperature ( K ).
1244 !~
1245 !~ References:
1246 !~ The algorithm is due to Nordquist, W.S., 1973: "Numerical approximations
1247 !~ of selected meteorological parameters for cloud physics problems,"
1248 !~ ecom-5475, Atmospheric Sciences Laboratory, U.S. Army Electronics
1249 !~ Command, White Sands Missile Range, New Mexico, 88002
1250 !~
1251 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1252 FUNCTION calc_es ( tK ) result ( es )
1254 implicit none
1256 !~ Variable Declaration
1257 ! --------------------
1258 real, intent ( in ) :: tK
1259 real :: es
1261 real :: p1, p2, c1
1263 p1 = 11.344 - 0.0303998 * tK
1264 p2 = 3.49149 - 1302.8844 / tK
1265 c1 = 23.832241 - 5.02808 * ALOG10 ( tK )
1266 es = 10.**(c1-1.3816E-7*10.**p1+8.1328E-3*10.**p2-2949.076/tK)
1268 END FUNCTION calc_es
1272 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1273 !~ ~!
1274 !~ Name: ~!
1275 !~ CATTurbulence ~!
1276 !~ ~!
1277 !~ Description: ~!
1278 !~ This function calculates the turbulence index ( TI ) for one layer ~!
1279 !~ in the atmosphere given the horizontal wind components and the geo- ~!
1280 !~ potential height of two pressure levels. The index is computed for ~!
1281 !~ the layer between the levels using the deformation and convergence ~!
1282 !~ of the wind field at the top and bottom of the layer and the vertical ~!
1283 !~ wind shear is calculated within the layer. The equation used for ~!
1284 !~ calculating TI is given by: ~!
1285 !~ ~!
1286 !~ ~!
1287 !~ TI = VWS * ( DEF + CONV ) ~!
1288 !~ ~!
1289 !~ Where: ~!
1290 !~ VWS = Vertical wind shear ~!
1291 !~ DEF = Deformation ~!
1292 !~ CONV = Convergence ~!
1293 !~ ~!
1294 !~ Notes: ~!
1295 !~ ~!
1296 !~ References: ~!
1297 !~ Ellrod, G.P. and D.J. Knapp, An objective clear-air turbulence ~!
1298 !~ forecasting technique: verification and operational use, Weather ~!
1299 !~ and Forecasting, 7, March 1992, pp. 150-165. ~!
1300 !~ ~!
1301 !~ Written: ~!
1302 !~ Scott Rentschler, Software Engineering Services ~!
1303 !~ Fine Scale Models Team ~!
1304 !~ Air Force Weather Agency, 16WS/WXN ~!
1305 !~ DSN: 271-3331 Comm: (402) 294-3331 ~!
1306 !~ scott.rentschler@offutt.af.mil ~!
1307 !~ ~!
1308 !~ History: ~!
1309 !~ 1 February 2008 ................... Scott Rentschler, (SES), 2WG/WEA ~!
1310 !~ INITIAL VERSION ~!
1311 !~ ~!
1312 !~ 8 July 2009 ....................... Scott Rentschler, (SES), 2WG/WEA ~!
1313 !~ Adapted for new driver. ~!
1314 !~ ~!
1315 !~ 1 November 2012 ......................... Scott Rentschler, 16WS/WXN ~!
1316 !~ Modified to accept layer argument, which adds the flexibility to make ~!
1317 !~ the computation for whichever flight level is desired. Cleaned up ~!
1318 !~ some of the code and added a couple comments. ~!
1319 !~ ~!
1320 !~ 28 August 2013 .................... Braedi Wickard, SEMS/NG/16WS/WXN ~!
1321 !~ Adapted for use within the Unified Post Processor. UPP can not handle ~!
1322 !~ the layer argument for flight levels, so reverted to hardcoded levels ~!
1323 !~ ~!
1324 !~ 25 April 2014 ............................. Glenn Creighton, 16WS/WXN ~!
1325 !~ Adapted for use within WRF. WRF already computes many of these terms. ~!
1326 !~ Stripped everything down to its bare bones to remove need to compute ~!
1327 !~ horizontal terms, now using deformation variables already within WRF. ~!
1328 !~ ~!
1329 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1330 FUNCTION CATTurbulence ( ugrdbot, ugrdtop, vgrdbot, vgrdtop &
1331 ,defor11bot, defor11top, defor12bot, defor12top &
1332 ,defor22bot, defor22top, zbot, ztop ) result ( ti )
1334 IMPLICIT NONE
1336 !~ Variable declarations
1337 ! ---------------------
1338 REAL, INTENT ( IN ) :: ugrdbot !~ U-wind bottom of layer
1339 REAL, INTENT ( IN ) :: ugrdtop !~ U-wind top of layer
1340 REAL, INTENT ( IN ) :: vgrdbot !~ V-wind bottom of layer
1341 REAL, INTENT ( IN ) :: vgrdtop !~ V-wind top of layer
1342 REAL, INTENT ( IN ) :: defor11bot !~ 2*du/dx at bottom of layer
1343 REAL, INTENT ( IN ) :: defor11top !~ 2*du/dx at top of layer
1344 REAL, INTENT ( IN ) :: defor12bot !~ du/dy+dv/dx at bottom of layer
1345 REAL, INTENT ( IN ) :: defor12top !~ du/dy+dv/dx at top of layer
1346 REAL, INTENT ( IN ) :: defor22bot !~ 2*dv/dy at bottom of layer
1347 REAL, INTENT ( IN ) :: defor22top !~ 2*dv/dy at top of layer
1348 REAL, INTENT ( IN ) :: zbot !~ Height grid bottom
1349 REAL, INTENT ( IN ) :: ztop !~ Height grid top
1350 REAL :: ti !~ Turbulence index
1352 !~ Function utility variables
1353 ! --------------------------
1354 REAL :: dudx, dudx1, dudx2 !~ Wind differentials
1355 REAL :: dvdy, dvdy1, dvdy2
1356 REAL :: dudz, dvdz
1358 REAL :: depth, vws, conv !~ Depth, vertical wind shear, convergence
1359 REAL :: def, shear, stretch !~ Deformation, shear, stretching terms
1361 !~ Initialize variables.
1362 ! ----------------------
1363 ti = REAL ( 0 )
1365 !~ Compute vertical wind shear
1366 ! ---------------------------
1367 depth = ABS ( ztop - zbot )
1368 dudz = ( ugrdbot - ugrdtop ) / depth
1369 dvdz = ( vgrdbot - vgrdtop ) / depth
1370 vws = SQRT ( dudz**2 + dvdz**2 )
1372 dudx1 = defor11top / 2.
1373 dudx2 = defor11bot / 2.
1374 dudx = ( dudx1 + dudx2 ) / REAL ( 2 )
1376 dvdy1 = defor22top / 2.
1377 dvdy2 = defor22bot / 2.
1378 dvdy = ( dvdy1 + dvdy2 ) / REAL ( 2 )
1380 !~ Compute the deformation
1381 ! -----------------------
1382 stretch = dudx - dvdy
1383 shear = ( defor12top + defor12bot ) / REAL ( 2 )
1384 def = SQRT ( stretch**2 + shear**2 )
1386 !~ Compute the convergence
1387 ! -----------------------
1388 conv = - ( dudx + dvdy )
1390 !~ Compute the turbulence index
1391 ! ----------------------------
1392 ti = vws * ( def + conv ) * 1.0E+07
1394 IF ( ti /= ti ) ti = REAL ( 0 )
1395 IF ( ti < 0 ) ti = REAL ( 0 )
1397 END FUNCTION CATTurbulence
1401 FUNCTION lin_interp ( x, f, y ) result ( g )
1403 ! Purpose:
1404 ! interpolates f(x) to point y
1405 ! assuming f(x) = f(x0) + a * (x - x0)
1406 ! where a = ( f(x1) - f(x0) ) / (x1 - x0)
1407 ! x0 <= x <= x1
1408 ! assumes x is monotonically increasing
1410 ! Author: D. Fillmore :: J. Done changed from r8 to r4
1411 ! Pilfered for AFWA diagnostics - G Creighton
1413 implicit none
1415 real, intent(in), dimension(:) :: x ! grid points
1416 real, intent(in), dimension(:) :: f ! grid function values
1417 real, intent(in) :: y ! interpolation point
1418 real :: g ! interpolated function value
1420 integer :: k ! interpolation point index
1421 integer :: n ! length of x
1422 real :: a
1424 n = size(x)
1426 ! find k such that x(k) < y =< x(k+1)
1427 ! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
1429 if (y <= x(1)) then
1430 k = 1
1431 else if (y >= x(n)) then
1432 k = n - 1
1433 else
1434 k = 1
1435 do while (y > x(k+1) .and. k < n)
1436 k = k + 1
1437 end do
1438 end if
1439 ! interpolate
1440 a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
1441 g = f(k) + a * (y - x(k))
1443 END FUNCTION lin_interp
1447 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1448 !~ ~!
1449 !~ Name: ~!
1450 !~ LLT_Windspeed ~!
1451 !~ ~!
1452 !~ Description: ~!
1453 !~ This function computes the dynamic term for the low-level turbulence ~!
1454 !~ algorithm. ~!
1455 !~ ~!
1456 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1457 FUNCTION LLT_Windspeed ( nlayer, u, v ) RESULT ( dynamic )
1458 IMPLICIT NONE
1460 !~ Variable Declaration
1461 ! --------------------
1462 INTEGER, INTENT ( IN ) :: nlayer
1463 REAL, INTENT ( IN ) :: u ( nlayer )
1464 REAL, INTENT ( IN ) :: v ( nlayer )
1465 REAL :: dynamic
1467 !~ Internal function variables
1468 ! ---------------------------
1469 INTEGER :: i
1470 REAL :: this_windspeed ( nlayer )
1471 REAL :: PI
1472 PARAMETER ( PI = 3.14159265359 )
1474 ! -------------------------------------------------------------------- !
1475 ! -------------------------------------------------------------------- !
1476 !~ Initialize variables
1477 ! --------------------
1478 dynamic = REAL ( 0 )
1480 !~ Compute the windspeed
1481 ! ---------------------
1482 DO i = 1, nlayer
1483 this_windspeed ( i ) = SQRT ( u(i)**2 + v(i)**2 )
1484 END DO
1486 !~ Compute the dynamic term
1487 ! -------------------------
1488 dynamic = ( this_windspeed(1)+this_windspeed(nlayer) ) / REAL (20)
1489 IF ( dynamic > REAL (2) ) dynamic = REAL ( 2 )
1490 dynamic = ( dynamic + REAL (3) ) / REAL ( 2 )
1491 dynamic = SIN ( dynamic*PI )
1492 dynamic = ( dynamic + REAL (1) ) / REAL ( 2 )
1495 END FUNCTION LLT_Windspeed
1498 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1499 !~ ~!
1500 !~ Name: ~!
1501 !~ LLT_Thermodynamic ~!
1502 !~ ~!
1503 !~ Description: ~!
1504 !~ This function computes the thermodynamic term for the low-level ~!
1505 !~ turbulence algorithm. ~!
1506 !~ ~!
1507 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1508 FUNCTION LLT_Thermodynamic ( nlayer, tK, hgt ) RESULT ( thermodynamic )
1509 IMPLICIT NONE
1511 !~ Variable Declaration
1512 ! --------------------
1513 INTEGER, INTENT ( IN ) :: nlayer
1514 REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperature (K)
1515 REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
1516 REAL :: thermodynamic
1518 !~ Internal function variables
1519 ! ---------------------------
1520 INTEGER :: i
1521 REAL :: lapse
1522 REAL :: PI
1523 PARAMETER ( PI = 3.14159265359 )
1525 ! -------------------------------------------------------------------- !
1526 ! -------------------------------------------------------------------- !
1528 !~ Initialize variables
1529 ! --------------------
1530 thermodynamic = REAL ( 0 )
1532 !~ Compute the lapse rate over the layer. The sign gets goofy here,
1533 !~ but works as coded below.
1534 ! -----------------------------------------------------------------
1535 lapse = ( tk(1) - tk(nlayer) ) * REAL ( 1000 )
1536 lapse = lapse / ( hgt(nlayer) - hgt(1) )
1538 !~ Compute the thermodynamic component
1539 ! -----------------------------------
1540 thermodynamic = lapse / REAL ( 10 )
1541 thermodynamic = ( thermodynamic + REAL (3) ) / REAL ( 2 )
1542 thermodynamic = SIN ( thermodynamic * PI )
1543 thermodynamic = ( thermodynamic + REAL (1) ) / REAL ( 2 )
1545 END FUNCTION LLT_Thermodynamic
1548 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1549 !~ ~!
1550 !~ Name: ~!
1551 !~ LLT_MountainWave ~!
1552 !~ ~!
1553 !~ Description: ~!
1554 !~ This function computes the mountain wave term for the low-level ~!
1555 !~ turbulence algorithm. ~!
1556 !~ ~!
1557 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1558 FUNCTION LLT_MountainWave ( nlayer, tdx, tdy, u, v, tK, hgt) &
1559 RESULT ( MountainWave )
1560 IMPLICIT NONE
1562 !~ Variable Declaration
1563 ! --------------------
1564 INTEGER, INTENT ( IN ) :: nlayer
1565 REAL, INTENT ( IN ) :: tdx !~ Terrain dx
1566 REAL, INTENT ( IN ) :: tdy !~ Terrain dy
1567 REAL, INTENT ( IN ) :: u ( nlayer ) !~ U components f
1568 REAL, INTENT ( IN ) :: v ( nlayer ) !~ V components
1569 REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperatures (K)
1570 REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
1571 REAL :: MountainWave !~ Mountain wave term
1573 !~ Internal function variables
1574 ! ---------------------------
1575 REAL :: u_term
1576 REAL :: v_term
1577 REAL :: uv_term
1578 REAL :: lapse
1579 REAL :: total_mw, this_total_mw
1580 REAL :: this_uv_term
1581 REAL :: min_uv_term, cross_terrain, max_total_mw
1582 INTEGER :: i, j, k
1584 REAL :: PI
1585 PARAMETER ( PI = 3.14159265359 )
1587 ! -------------------------------------------------------------------- !
1588 ! -------------------------------------------------------------------- !
1590 !~ Initialize variables
1591 ! --------------------
1592 MountainWave = REAL ( 0 )
1594 !~ Loop through the layer
1595 ! ----------------------
1596 DO i = 2, nlayer
1598 !~ Wind terrain term
1599 ! -----------------
1600 u_term = ( (u(i-1) + u(i) ) / REAL(2) ) * tdx
1601 v_term = ( (v(i-1) + v(i) ) / REAL(2) ) * tdy
1602 this_uv_term = ( u_term + v_term ) * REAL ( -1 )
1603 !IF ( uv_term < REAL (0) ) uv_term = REAL ( 0 )
1604 IF ( min_uv_term < REAL (0) ) min_uv_term = REAL ( 0 )
1605 IF ( i == 2 ) THEN
1606 !uv_term = this_uv_term
1607 min_uv_term = this_uv_term
1608 ELSE
1609 !IF ( this_uv_term < uv_term ) uv_term = this_uv_term
1610 IF ( this_uv_term < min_uv_term ) min_uv_term = this_uv_term
1611 END IF
1613 !~ Lapse rate
1614 ! ----------
1615 lapse = ( tK (i-1) - tK (i) ) * REAL ( 1000 )
1616 lapse = lapse / ABS ( hgt(i)-hgt(i-1) )
1617 IF ( lapse > REAL (0) ) lapse = REAL ( 0 )
1618 lapse = lapse * REAL ( -1 )
1620 this_total_mw = this_uv_term * lapse * REAL ( 40000 )
1621 IF ( i == 2 ) THEN
1622 total_mw = this_total_mw
1623 ELSE
1624 IF ( this_total_mw > total_mw ) total_mw = this_total_mw
1625 END IF
1627 END DO
1629 !min_uv_term = uv_term
1630 cross_terrain = min_uv_term * REAL ( 500 )
1632 IF ( min_uv_term < 0.03 ) THEN
1633 cross_terrain = REAL ( 0 )
1634 END IF
1636 IF ( cross_terrain > REAL (50) ) cross_terrain = REAL ( 50 )
1638 !~ Multiply the lapse (inversion) array and the mountain wave array
1639 ! ----------------------------------------------------------------
1640 IF ( total_mw > REAL (50) ) total_mw = REAL ( 50 )
1642 !~ Add the cross terrain flow and inversion term
1643 ! ---------------------------------------------
1644 MountainWave = ( total_mw*(cross_terrain/50.) ) + cross_terrain
1645 MountainWave = MountainWave / REAL ( 100 )
1647 END FUNCTION LLT_MountainWave
1651 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1652 !~ ~!
1653 !~ Name: ~!
1654 !~ LLT_TrappedWave ~!
1655 !~ ~!
1656 !~ Description: ~!
1657 !~ This function computes the trapped wave term for the low-level ~!
1658 !~ turbulence algorithm. ~!
1659 !~ ~!
1660 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1661 FUNCTION LLT_TrappedWave ( nlayer, u, v, p ) RESULT ( trapped )
1662 IMPLICIT NONE
1664 !~ Variable Declaration
1665 ! --------------------
1666 INTEGER, INTENT ( IN ) :: nlayer
1667 REAL, INTENT ( IN ) :: u ( nlayer )
1668 REAL, INTENT ( IN ) :: v ( nlayer )
1669 REAL, INTENT ( IN ) :: p ( nlayer )
1670 REAL :: trapped
1672 !~ Internal function variables
1673 ! ---------------------------
1674 INTEGER :: i
1675 REAL :: du, dv
1676 REAL :: scale_fact, this_p
1677 REAL :: dudv, this_dudv
1678 REAL :: PI
1679 PARAMETER ( PI = 3.14159265359 )
1681 !~ Scale parameters
1682 ! ----------------
1683 REAL, PARAMETER :: scale_950 = 0.050000 !~ 1/20
1684 REAL, PARAMETER :: scale_925 = 0.040000 !~ 1/25
1685 REAL, PARAMETER :: scale_900 = 0.025000 !~ 1/40
1686 REAL, PARAMETER :: scale_850 = 0.010000 !~ 1/100
1687 REAL, PARAMETER :: scale_800 = 0.005000 !~ 1/200
1688 REAL, PARAMETER :: scale_750 = 0.002941 !~ 1/340
1689 REAL, PARAMETER :: scale_700 = 0.001923 !~ 1/520
1690 REAL, PARAMETER :: scale_650 = 0.001351 !~ 1/740
1691 REAL, PARAMETER :: scale_600 = 0.001000 !~ 1/1000
1692 REAL, PARAMETER :: scale_550 = 0.000800 !~ 1/1250
1694 ! -------------------------------------------------------------------- !
1695 ! -------------------------------------------------------------------- !
1697 !~ Initialize variables
1698 ! --------------------
1699 trapped = REAL ( 0 )
1701 !~ Compute the trapped wave term
1702 ! ------------------
1703 dudv = REAL ( 0 )
1704 DO i = 2, nlayer
1706 !~ Compute dudv first
1707 ! ------------------
1708 du = u ( i-1 ) - u ( i )
1709 dv = v ( i-1 ) - v ( i )
1711 !~ Scale based on pressure level
1712 ! -----------------------------
1713 this_p = p ( i ) / REAL ( 100 )
1714 IF ( this_p > REAL (950) ) THEN
1715 scale_fact = scale_950
1716 ELSE IF ( this_p <= REAL (950) .AND. this_p > REAL (925) ) THEN
1717 scale_fact = scale_925
1718 ELSE IF ( this_p <= REAL (925) .AND. this_p > REAL (900) ) THEN
1719 scale_fact = scale_900
1720 ELSE IF ( this_p <= REAL (900) .AND. this_p > REAL (850) ) THEN
1721 scale_fact = scale_850
1722 ELSE IF ( this_p <= REAL (850) .AND. this_p > REAL (800) ) THEN
1723 scale_fact = scale_800
1724 ELSE IF ( this_p <= REAL (800) .AND. this_p > REAL (750) ) THEN
1725 scale_fact = scale_750
1726 ELSE IF ( this_p <= REAL (750) .AND. this_p > REAL (700) ) THEN
1727 scale_fact = scale_700
1728 ELSE IF ( this_p <= REAL (700) .AND. this_p > REAL (650) ) THEN
1729 scale_fact = scale_650
1730 ELSE IF ( this_p <= REAL (650) .AND. this_p > REAL (600) ) THEN
1731 scale_fact = scale_600
1732 ELSE IF ( this_p <= REAL (600) ) THEN
1733 scale_fact = scale_550
1734 END IF
1736 this_dudv = ( (du**2)*(dv**2) ) * scale_fact
1737 IF ( this_dudv > dudv ) dudv = this_dudv
1739 END DO
1741 trapped = dudv
1742 IF ( trapped > REAL ( 1 ) ) trapped = REAL ( 1 )
1743 trapped = trapped / REAL ( 4 )
1745 END FUNCTION LLT_TrappedWave
1747 END MODULE diag_functions
1748 #endif