| blob | 0466dff45f90c6e817bfbfc5ceec34dc97f142ef |
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) .and. (lfclev >= 0)) THEN !Modified by Zhixiao
764 plfc = p (lfclev)
765 pel = p (ellev)
766 CIN = REAL ( 0 )
767 DO k = sfc+1, nz
768 ! Make CAPE and CIN consistent with AMS definition
769 ! https://glossary.ametsoc.org/wiki/Convective_available_potential_energy
770 ! https://glossary.ametsoc.org/wiki/Convective_inhibition
771 IF ( p (k) <= plcl .and. p (k) > plfc) THEN !Modified by Zhixiao
772 ! CIN is the vertically integrated negative buoyant energy between LCL and LFC
773 CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
774 END IF
775 IF ( p (k) <= plfc .and. p (k) > pel) THEN !Modified by Zhixiao
776 ! CAPE is the vertically integrated positive buoyant energy between LFC and EL
777 CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
778 END IF !Modified by Zhixiao
779 END DO !Modified by Zhixiao
780 END IF !Modified by Zhixiao
781 !~ Calculate lifted index by interpolating difference between
782 !~ parcel and ambient Tv to 500mb.
783 ! ----------------------------------------------------------
784 DO k = sfc + 1, nz
786 pm = 50000.
787 pu = p ( k )
788 pd = p ( k - 1 )
790 !~ If we're already above 500mb just set lifted index to 0.
791 !~ --------------------------------------------------------
792 IF ( pd .le. pm ) THEN
793 lidx = 0.
794 EXIT
796 ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
798 !~ Found trapping pressure: up, middle, down.
799 !~ We are doing first order interpolation.
800 ! ------------------------------------------
801 lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
802 lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
803 lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
804 EXIT
806 ENDIF
808 END DO
809 !~ Assuming the the LFC is at a pressure level for now
810 ! ---------------------------------------------------
811 IF ( lfclev > 0 ) THEN
812 PLFC = p ( lfclev )
813 ZLFC = hgt ( lfclev )
814 END IF
816 IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN
817 PLFC = REAL ( -1 )
818 ZLFC = REAL ( -1 )
819 END IF
821 IF ( CAPE /= CAPE ) cape = REAL ( 0 )
823 IF ( CIN /= CIN ) cin = RUNDEF
825 !~ Chirp
826 ! -----
827 ! WRITE ( *,* ) ' CAPE: ', cape, ' CIN: ', cin
828 ! WRITE ( *,* ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
829 ! WRITE ( *,* ) ''
830 ! WRITE ( *,* ) ' Exiting buoyancy.'
831 ! WRITE ( *,* ) ' ==================================== '
832 ! WRITE ( *,* ) ''
834 END FUNCTION Buoyancy
838 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
839 ! . . .
840 ! SUBPROGRAM: NGMSLP NMC SEA LEVEL PRESSURE REDUCTION
841 ! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-02
842 !
843 ! ABSTRACT:
844 !
845 ! THIS ROUTINE COMPUTES SEA LEVEL PRESSURE USING THE
846 ! HYDROSTATIC EQUATION WITH THE SHUELL CORRECTION. THE
847 ! FOLLOWING IS BASED ON DOCUMENTATION IN SUBROUTINE
848 ! OUTHYDRO OF THE NGM:
849 !
850 ! THE FUNDAMENTAL HYDROSTATIC EQUATION IS
851 ! D(HEIGHT)
852 ! --------- = TAU = VIRTUAL TEMPERATURE * (RGAS/GRAVITY)
853 ! D (Z)
854 ! WHERE
855 ! Z = MINUS LOG OF PRESSURE (-LN(P)).
856 !
857 ! SEA-LEVEL PRESSURE IS COMPUTED FROM THE FORMULA
858 ! PRESS(MSL) = PRESS(GROUND) * EXP( F)
859 ! WHERE
860 ! F = HEIGHT OF GROUND / MEAN TAU
861 ! MEAN TAU = ( TAU(GRND) + TAU(SL) ) / 2
862 !
863 ! IN THE NGM TAU(GRND) AND TAU(SL) ARE FIRST SET USING A
864 ! 6.5DEG/KM LAPSE RATE FROM LOWEST MDL LEVEL. THIS IS MODIFIED
865 ! BY A CORRECTION BASED ON THE CRITICAL TAU OF THE SHUELL
866 ! CORRECTION:
867 ! TAUCR=(RGASD/GRAVITY) * 290.66
868 !
869 ! 1) WHERE ONLY TAU(SL) EXCEEDS TAUCR, CHANGE TAU(SL) TO TAUCR.
870 !
871 ! 2) WHERE BOTH TAU(SL) AND TAU(GRND) EXCEED TAUCR,
872 ! CHANGE TAU(SL) TO TAUCR-CONST*(TAU(GRND)-TAUCR )**2
873 ! WHERE CONST = .005 (GRAVITY/RGASD)
874 !
875 ! THE AVERAGE OF TAU(SL) AND TAU(GRND) IS THEN USED TOGETHER
876 ! WITH THE GROUND HEIGHT AND PRESSURE TO DERIVE THE PRESSURE
877 ! AT SEA LEVEL.
878 !
879 ! HEIGHT OF THE 1000MB SURFACE IS COMPUTED FROM THE MSL PRESSURE
880 ! FIELD USING THE FORMULA:
881 !
882 ! P(MSL) - P(1000MB) = MEAN DENSITY * GRAVITY * HGT(1000MBS)
883 !
884 ! WHERE P(MSL) IS THE SEA LEVEL PRESSURE FIELD WE HAVE JUST
885 ! COMPUTED.
886 !
887 !
888 ! MEB 6/13/02: THIS CODE HAS BEEN SIMPLIFIED CONSIDERABLY FROM
889 ! THE ONE USED IN ETAPOST. HORIZONTAL SMOOTHING HAS BEEN
890 ! REMOVED AND THE FIRST MODEL LEVEL IS USED RATHER
891 ! THAN THE MEAN OF THE VIRTUAL TEMPERATURES IN
892 ! THE LOWEST 30MB ABOVE GROUND TO COMPUTE TAU(GRND).
893 !
894 ! .
895 !
896 ! PROGRAM HISTORY LOG:
897 ! 93-02-02 RUSS TREADON
898 ! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D
899 ! 00-01-04 JIM TUCCILLO - MPI VERSION
900 ! 01-10-25 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
901 ! 01-11-02 H CHUANG - MODIFIED LINE 234 FOR COMPUTATION OF
902 ! SIGMA/HYBRID SLP
903 ! 01-12-18 H CHUANG - INCLUDED SMOOTHING ALONG BOUNDARIES TO BE
904 ! CONSISTENT WITH MESINGER SLP
905 ! 02-06-13 MIKE BALDWIN - WRF VERSION
906 ! 06-12-18 H CHUANG - BUG FIX TO CORRECT TAU AT SFC
907 ! 14-04-17 G CREIGHTON - MODIFIED TO INSERT INTO AFWA DIAGNOSTICS IN WRF
908 !
909 !$$$
911 FUNCTION MSLP ( zsfc, psfc, zlev1, qlev1, tlev1 )
913 implicit none
916 ! DECLARE VARIABLES
918 REAL, INTENT ( IN ) :: zsfc !~ Surface height ( m )
919 REAL, INTENT ( IN ) :: psfc !~ Surface height ( m )
920 REAL, INTENT ( IN ) :: zlev1 !~ Level 1 height ( m )
921 REAL, INTENT ( IN ) :: qlev1 !~ Level 1 mixing ratio ( kg/kg )
922 REAL, INTENT ( IN ) :: tlev1 !~ Level 1 temperature ( K )
923 real,PARAMETER :: G=9.81
924 real,PARAMETER :: GI=1./G
925 real,PARAMETER :: RD=287.0
926 real,PARAMETER :: ZSL=0.0
927 real,PARAMETER :: TAUCR=RD*GI*290.66,CONST=0.005*G/RD
928 real,PARAMETER :: GORD=G/RD,DP=60.E2
929 real,PARAMETER :: GAMMA=6.5E-3
931 real MSLP,TVRT,TVRSFC,TAUSFC,TVRSL,TAUSL,TAUAVG
932 !
933 !**********************************************************************
934 ! START NGMSLP HERE.
935 !
936 MSLP = PSFC
937 !
938 ! COMPUTE LAYER TAU (VIRTUAL TEMP*RD/G).
939 TVRT = TLEV1*(1.0+0.608*QLEV1)
940 !TAU = TVRT*RD*GI
941 !
942 ! COMPUTE TAU AT THE GROUND (Z=ZSFC) AND SEA LEVEL (Z=0)
943 ! ASSUMING A CONSTANT LAPSE RATE OF GAMMA=6.5DEG/KM.
944 TVRSFC = TVRT + (ZLEV1 - ZSFC)*GAMMA
945 TAUSFC = TVRSFC*RD*GI
946 TVRSL = TVRT + (ZLEV1 - ZSL)*GAMMA
947 TAUSL = TVRSL*RD*GI
948 !
949 ! IF NEED BE APPLY SHEULL CORRECTION.
950 IF ((TAUSL.GT.TAUCR).AND.(TAUSFC.LE.TAUCR)) THEN
951 TAUSL=TAUCR
952 ELSEIF ((TAUSL.GT.TAUCR).AND.(TAUSFC.GT.TAUCR)) THEN
953 TAUSL = TAUCR-CONST*(TAUSFC-TAUCR)**2
954 ENDIF
955 !
956 ! COMPUTE MEAN TAU.
957 TAUAVG = 0.5*(TAUSL+TAUSFC)
958 !
959 ! COMPUTE SEA LEVEL PRESSURE.
960 MSLP = PSFC*EXP(ZSFC/TAUAVG)
962 END FUNCTION MSLP
966 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
967 !~ ~!
968 !~ Name: ~!
969 !~ calc_fits ~!
970 !~ ~!
971 !~ Description: ~!
972 !~ This function computes Fighter Index Thermal Stress values given ~!
973 !~ dry bulb temperature, relative humidity, and pressure. ~!
974 !~ ~!
975 !~ Usage: ~!
976 !~ fitsval = calc_fits ( p, tK, rh ) ~!
977 !~ ~!
978 !~ Where: ~!
979 !~ p = Pressure ( Pa ) ~!
980 !~ tK = Temperature ( K ) ~!
981 !~ rh = Relative Humidity ( % ) ~!
982 !~ ~!
983 !~ Reference: ~!
984 !~ Stribley, R.F., S. Nunneley, 1978: Fighter Index of Thermal Stress: ~!
985 !~ Development of interim guidance for hot-weather USAF operations. ~!
986 !~ SAM-TR-78-6. Eqn. 9 ~!
987 !~ ~!
988 !~ Formula: ~!
989 !~ FITS = 0.8281*Twb + 0.3549*Tdb + 5.08 (degrees Celsius) ~!
990 !~ ~!
991 !~ Where: ~!
992 !~ Twb = Wet Bulb Temperature ~!
993 !~ Tdb = Dry Bulb Temperature ~!
994 !~ ~!
995 !~ Written: ~!
996 !~ Scott Rentschler, Software Engineering Services ~!
997 !~ Fine Scale Models Team ~!
998 !~ Air Force Weather Agency, 16WS/WXN ~!
999 !~ DSN: 271-3331 Comm: (402) 294-3331 ~!
1000 !~ scott.rentschler@offutt.af.mil ~!
1001 !~ ~!
1002 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1003 FUNCTION calc_fits ( p, tK, rh ) RESULT ( fits )
1005 implicit none
1007 !~ Variable declaration
1008 ! --------------------
1009 real, intent ( in ) :: p !~ Pressure ( Pa )
1010 real, intent ( in ) :: tK !~ Temperature ( K )
1011 real, intent ( in ) :: rh !~ Rel Humidity ( % )
1012 real :: fits !~ FITS index value
1014 !~ Utility variables
1015 ! --------------------------
1016 real :: twb !~ Wet bulb temperature ( C )
1017 real :: wbt
1019 ! ---------------------------------------------------------------------- !
1020 ! ---------------------------------------------------------------------- !
1022 !~ Initialize variables
1023 ! --------------------
1024 fits = REAL ( 0 )
1026 !~ Get the wet bulb temperature in degrees Celsius
1027 ! -----------------------------------------------
1028 twb = WetBulbTemp ( p, tK, rh ) - 273.15
1030 !~ Compute the FITS index
1031 ! ----------------------
1032 fits = 0.8281*twb + 0.3549*( tK - 273.15 ) + 5.08
1034 !~ Convert the index to Kelvin
1035 ! ---------------------------
1036 fits = fits + 273.15
1038 END FUNCTION calc_fits
1042 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1043 !~
1044 !~ Name:
1045 !~ calc_wc
1046 !~
1047 !~ Description:
1048 !~ This function calculates wind chill given temperature ( K ) and
1049 !~ wind speed ( m/s )
1050 !~
1051 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1052 FUNCTION calc_wc ( tK, wspd ) RESULT ( wc )
1054 implicit none
1056 !~ Variable Declarations
1057 ! ---------------------
1058 real, intent ( in ) :: tK
1059 real, intent ( in ) :: wspd
1061 real :: tF, wc, wspd_mph
1063 wspd_mph = wspd * 2.23693629 ! convert to mph
1064 tF = (( tK - 273.15 ) * ( REAL (9) / REAL (5) ) ) + REAL ( 32 )
1066 wc = 35.74 &
1067 + ( 0.6215 * tF ) &
1068 - ( 35.75 * ( wspd_mph**0.16 ) ) &
1069 + ( 0.4275 * tF * ( wspd_mph**0.16 ) )
1071 wc = (( wc - REAL (32) ) * ( REAL (5) / REAL (9) ) ) + 273.15
1073 END FUNCTION calc_wc
1077 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1078 !~
1079 !~ Name:
1080 !~ calc_hi
1081 !~
1082 !~ Description:
1083 !~ This subroutine calculates the heat index. Requires temperature ( K )
1084 !~ and relative humidity ( % ).
1085 !~
1086 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1087 FUNCTION calc_hi ( Tk, RH ) result ( HI )
1089 implicit none
1091 !~ Variable declarations
1092 ! ---------------------
1093 real, intent ( in ) :: Tk
1094 real, intent ( in ) :: RH
1096 real :: tF, tF2, rh2, HI
1098 !~ If temperature > 70F then calculate heat index, else set it equal
1099 !~ to dry temperature
1100 ! -----------------------------------------------------------------
1101 IF ( Tk > 294.26111 ) THEN
1103 tF = ( (Tk - 273.15) * (REAL (9)/REAL (5)) ) + REAL ( 32 )
1104 tF2 = tF ** 2
1105 rh2 = RH ** 2
1107 HI = -42.379 &
1108 + ( 2.04901523 * tF ) &
1109 + ( 10.14333127 * RH ) &
1110 - ( 0.22475541 * tF * RH ) &
1111 - ( 6.83783E-03 * tF2 ) &
1112 - ( 5.481717E-02 * rh2 ) &
1113 + ( 1.22874E-03 * tF2 * RH ) &
1114 + ( 8.5282E-04 * tF * rh2 ) &
1115 - ( 1.99E-06 * tF2 * rh2 )
1117 HI = ((HI - REAL (32)) * (REAL (5)/REAL (9))) + 273.15
1118 ELSE
1119 HI = Tk
1120 END IF
1122 END FUNCTION calc_hi
1124 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1125 !~ ~!
1126 !~ Name: ~!
1127 !~ WetBulbTemp ~!
1128 !~ ~!
1129 !~ Description: ~!
1130 !~ This function approximates the Wet Bulb Temperature (K) provided ~!
1131 !~ dry bulb temperature (K), relative humidity (%), and pressure (Pa). ~!
1132 !~ ~!
1133 !~ Usage: ~!
1134 !~ wbt = WetBulbTemperature ( p, tK, rh ) ~!
1135 !~ ~!
1136 !~ Where: ~!
1137 !~ p = Pressure ( Pa ) ~!
1138 !~ tK = Temperature ( K ) ~!
1139 !~ rh = Relative Humidity ( % ) ~!
1140 !~ ~!
1141 !~ Reference: ~!
1142 !~ American Society of Civil Engineers ~!
1143 !~ Evapotraspiration and Irrigation Water Requirements ~!
1144 !~ Jensen et al (1990) ASCE Manual No. 70, pp 176-177 ~!
1145 !~ ~!
1146 !~ Written: ~!
1147 !~ Scott Rentschler, Software Engineering Services ~!
1148 !~ Fine Scale Models Team ~!
1149 !~ Air Force Weather Agency ~!
1150 !~ DSM: 271-3331 Comm: (402) 294-3331 ~!
1151 !~ scott.rentschler@offutt.af.mil ~!
1152 !~ ~!
1153 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1154 FUNCTION WetBulbTemp ( p, tK, rh) result( wbt )
1156 implicit none
1158 !~ Variable delclaration
1159 ! ---------------------
1160 real, intent ( in ) :: p !~ Pressure ( Pa )
1161 real, intent ( in ) :: tK !~ Temperature ( K )
1162 real, intent ( in ) :: rh !~ Relative Humidity ( % )
1163 real :: wbt !~ Wet Bulb Temperature ( K )
1165 !~ Utility variables
1166 ! -----------------
1167 real :: tdK !~ Dewpoint temperature ( K )
1168 real :: tC !~ Temperature ( C )
1169 real :: tdC !~ Dewpoint temperature ( K )
1170 real :: svapr !~ Saturation vapor pressure ( Pa )
1171 real :: vapr !~ Ambient vapor pressure ( Pa )
1172 real :: gamma !~ Dummy term
1173 real :: delta !~ Dummy term
1175 ! ---------------------------------------------------------------------- !
1176 ! ---------------------------------------------------------------------- !
1177 !~ Initialize variables
1178 ! --------------------
1179 wbt = REAL ( 0 )
1180 tC = tK - 273.15
1182 !~ Compute saturation vapor pressure ( Pa )
1183 ! ----------------------------------------
1184 svapr = calc_es ( tK ) * REAL ( 100 )
1186 !~ Compute vapor pressure
1187 ! ----------------------
1188 vapr = svapr * ( rh / REAL (100) )
1190 !~ Grab the dewpoint
1191 ! -----------------
1192 tdC = calc_Dewpoint ( tC, rh )
1193 tdK = tdC + 273.15
1195 !~ Compute dummy terms
1196 ! -------------------
1197 gamma = 0.00066 * ( p / REAL (1000) )
1198 delta = REAL ( 4098 ) * ( vapr / REAL(1000) ) / ( (tC+237.3)**2 )
1200 !~ Compute the wet bulb temperature
1201 ! --------------------------------
1202 wbt = ( ((gamma * tC) + (delta * tdC)) / (gamma + delta) ) + 273.15
1204 END FUNCTION WetBulbTemp
1207 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1208 !~
1209 !~ Name:
1210 !~ calc_Dewpoint
1211 !~
1212 !~ Description:
1213 !~ This function approximates dewpoint given temperature and rh.
1214 !~
1215 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1216 FUNCTION calc_Dewpoint ( tC, rh) result( Dewpoint )
1218 implicit none
1220 !~ Variable Declaration
1221 ! --------------------
1222 real, intent ( in ) :: tC
1223 real, intent ( in ) :: rh
1224 real :: Dewpoint
1226 real :: term, es, e1, e, logs, expon
1228 expon = ( 7.5*tC ) / ( 237.7+tC )
1229 es = 6.112 * ( 10**expon ) ! Saturated vapor pressure
1230 e = es * ( rh/100.0 ) ! Vapor pressure
1231 logs = LOG10 ( e/6.112 )
1232 Dewpoint = ( 237.7*logs ) / ( 7.5-logs )
1234 END FUNCTION calc_Dewpoint
1237 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1238 !~
1239 !~ Name:
1240 !~ calc_es
1241 !~
1242 !~ Description:
1243 !~ This function returns the saturation vapor pressure over water ( hPa )
1244 !~ given temperature ( K ).
1245 !~
1246 !~ References:
1247 !~ The algorithm is due to Nordquist, W.S., 1973: "Numerical approximations
1248 !~ of selected meteorological parameters for cloud physics problems,"
1249 !~ ecom-5475, Atmospheric Sciences Laboratory, U.S. Army Electronics
1250 !~ Command, White Sands Missile Range, New Mexico, 88002
1251 !~
1252 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1253 FUNCTION calc_es ( tK ) result ( es )
1255 implicit none
1257 !~ Variable Declaration
1258 ! --------------------
1259 real, intent ( in ) :: tK
1260 real :: es
1262 real :: p1, p2, c1
1264 p1 = 11.344 - 0.0303998 * tK
1265 p2 = 3.49149 - 1302.8844 / tK
1266 c1 = 23.832241 - 5.02808 * ALOG10 ( tK )
1267 es = 10.**(c1-1.3816E-7*10.**p1+8.1328E-3*10.**p2-2949.076/tK)
1269 END FUNCTION calc_es
1273 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1274 !~ ~!
1275 !~ Name: ~!
1276 !~ CATTurbulence ~!
1277 !~ ~!
1278 !~ Description: ~!
1279 !~ This function calculates the turbulence index ( TI ) for one layer ~!
1280 !~ in the atmosphere given the horizontal wind components and the geo- ~!
1281 !~ potential height of two pressure levels. The index is computed for ~!
1282 !~ the layer between the levels using the deformation and convergence ~!
1283 !~ of the wind field at the top and bottom of the layer and the vertical ~!
1284 !~ wind shear is calculated within the layer. The equation used for ~!
1285 !~ calculating TI is given by: ~!
1286 !~ ~!
1287 !~ ~!
1288 !~ TI = VWS * ( DEF + CONV ) ~!
1289 !~ ~!
1290 !~ Where: ~!
1291 !~ VWS = Vertical wind shear ~!
1292 !~ DEF = Deformation ~!
1293 !~ CONV = Convergence ~!
1294 !~ ~!
1295 !~ Notes: ~!
1296 !~ ~!
1297 !~ References: ~!
1298 !~ Ellrod, G.P. and D.J. Knapp, An objective clear-air turbulence ~!
1299 !~ forecasting technique: verification and operational use, Weather ~!
1300 !~ and Forecasting, 7, March 1992, pp. 150-165. ~!
1301 !~ ~!
1302 !~ Written: ~!
1303 !~ Scott Rentschler, Software Engineering Services ~!
1304 !~ Fine Scale Models Team ~!
1305 !~ Air Force Weather Agency, 16WS/WXN ~!
1306 !~ DSN: 271-3331 Comm: (402) 294-3331 ~!
1307 !~ scott.rentschler@offutt.af.mil ~!
1308 !~ ~!
1309 !~ History: ~!
1310 !~ 1 February 2008 ................... Scott Rentschler, (SES), 2WG/WEA ~!
1311 !~ INITIAL VERSION ~!
1312 !~ ~!
1313 !~ 8 July 2009 ....................... Scott Rentschler, (SES), 2WG/WEA ~!
1314 !~ Adapted for new driver. ~!
1315 !~ ~!
1316 !~ 1 November 2012 ......................... Scott Rentschler, 16WS/WXN ~!
1317 !~ Modified to accept layer argument, which adds the flexibility to make ~!
1318 !~ the computation for whichever flight level is desired. Cleaned up ~!
1319 !~ some of the code and added a couple comments. ~!
1320 !~ ~!
1321 !~ 28 August 2013 .................... Braedi Wickard, SEMS/NG/16WS/WXN ~!
1322 !~ Adapted for use within the Unified Post Processor. UPP can not handle ~!
1323 !~ the layer argument for flight levels, so reverted to hardcoded levels ~!
1324 !~ ~!
1325 !~ 25 April 2014 ............................. Glenn Creighton, 16WS/WXN ~!
1326 !~ Adapted for use within WRF. WRF already computes many of these terms. ~!
1327 !~ Stripped everything down to its bare bones to remove need to compute ~!
1328 !~ horizontal terms, now using deformation variables already within WRF. ~!
1329 !~ ~!
1330 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1331 FUNCTION CATTurbulence ( ugrdbot, ugrdtop, vgrdbot, vgrdtop &
1332 ,defor11bot, defor11top, defor12bot, defor12top &
1333 ,defor22bot, defor22top, zbot, ztop ) result ( ti )
1335 IMPLICIT NONE
1337 !~ Variable declarations
1338 ! ---------------------
1339 REAL, INTENT ( IN ) :: ugrdbot !~ U-wind bottom of layer
1340 REAL, INTENT ( IN ) :: ugrdtop !~ U-wind top of layer
1341 REAL, INTENT ( IN ) :: vgrdbot !~ V-wind bottom of layer
1342 REAL, INTENT ( IN ) :: vgrdtop !~ V-wind top of layer
1343 REAL, INTENT ( IN ) :: defor11bot !~ 2*du/dx at bottom of layer
1344 REAL, INTENT ( IN ) :: defor11top !~ 2*du/dx at top of layer
1345 REAL, INTENT ( IN ) :: defor12bot !~ du/dy+dv/dx at bottom of layer
1346 REAL, INTENT ( IN ) :: defor12top !~ du/dy+dv/dx at top of layer
1347 REAL, INTENT ( IN ) :: defor22bot !~ 2*dv/dy at bottom of layer
1348 REAL, INTENT ( IN ) :: defor22top !~ 2*dv/dy at top of layer
1349 REAL, INTENT ( IN ) :: zbot !~ Height grid bottom
1350 REAL, INTENT ( IN ) :: ztop !~ Height grid top
1351 REAL :: ti !~ Turbulence index
1353 !~ Function utility variables
1354 ! --------------------------
1355 REAL :: dudx, dudx1, dudx2 !~ Wind differentials
1356 REAL :: dvdy, dvdy1, dvdy2
1357 REAL :: dudz, dvdz
1359 REAL :: depth, vws, conv !~ Depth, vertical wind shear, convergence
1360 REAL :: def, shear, stretch !~ Deformation, shear, stretching terms
1362 !~ Initialize variables.
1363 ! ----------------------
1364 ti = REAL ( 0 )
1366 !~ Compute vertical wind shear
1367 ! ---------------------------
1368 depth = ABS ( ztop - zbot )
1369 dudz = ( ugrdbot - ugrdtop ) / depth
1370 dvdz = ( vgrdbot - vgrdtop ) / depth
1371 vws = SQRT ( dudz**2 + dvdz**2 )
1373 dudx1 = defor11top / 2.
1374 dudx2 = defor11bot / 2.
1375 dudx = ( dudx1 + dudx2 ) / REAL ( 2 )
1377 dvdy1 = defor22top / 2.
1378 dvdy2 = defor22bot / 2.
1379 dvdy = ( dvdy1 + dvdy2 ) / REAL ( 2 )
1381 !~ Compute the deformation
1382 ! -----------------------
1383 stretch = dudx - dvdy
1384 shear = ( defor12top + defor12bot ) / REAL ( 2 )
1385 def = SQRT ( stretch**2 + shear**2 )
1387 !~ Compute the convergence
1388 ! -----------------------
1389 conv = - ( dudx + dvdy )
1391 !~ Compute the turbulence index
1392 ! ----------------------------
1393 ti = vws * ( def + conv ) * 1.0E+07
1395 IF ( ti /= ti ) ti = REAL ( 0 )
1396 IF ( ti < 0 ) ti = REAL ( 0 )
1398 END FUNCTION CATTurbulence
1402 FUNCTION lin_interp ( x, f, y ) result ( g )
1404 ! Purpose:
1405 ! interpolates f(x) to point y
1406 ! assuming f(x) = f(x0) + a * (x - x0)
1407 ! where a = ( f(x1) - f(x0) ) / (x1 - x0)
1408 ! x0 <= x <= x1
1409 ! assumes x is monotonically increasing
1411 ! Author: D. Fillmore :: J. Done changed from r8 to r4
1412 ! Pilfered for AFWA diagnostics - G Creighton
1414 implicit none
1416 real, intent(in), dimension(:) :: x ! grid points
1417 real, intent(in), dimension(:) :: f ! grid function values
1418 real, intent(in) :: y ! interpolation point
1419 real :: g ! interpolated function value
1421 integer :: k ! interpolation point index
1422 integer :: n ! length of x
1423 real :: a
1425 n = size(x)
1427 ! find k such that x(k) < y =< x(k+1)
1428 ! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
1430 if (y <= x(1)) then
1431 k = 1
1432 else if (y >= x(n)) then
1433 k = n - 1
1434 else
1435 k = 1
1436 do while (y > x(k+1) .and. k < n)
1437 k = k + 1
1438 end do
1439 end if
1440 ! interpolate
1441 a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
1442 g = f(k) + a * (y - x(k))
1444 END FUNCTION lin_interp
1448 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1449 !~ ~!
1450 !~ Name: ~!
1451 !~ LLT_Windspeed ~!
1452 !~ ~!
1453 !~ Description: ~!
1454 !~ This function computes the dynamic term for the low-level turbulence ~!
1455 !~ algorithm. ~!
1456 !~ ~!
1457 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1458 FUNCTION LLT_Windspeed ( nlayer, u, v ) RESULT ( dynamic )
1459 IMPLICIT NONE
1461 !~ Variable Declaration
1462 ! --------------------
1463 INTEGER, INTENT ( IN ) :: nlayer
1464 REAL, INTENT ( IN ) :: u ( nlayer )
1465 REAL, INTENT ( IN ) :: v ( nlayer )
1466 REAL :: dynamic
1468 !~ Internal function variables
1469 ! ---------------------------
1470 INTEGER :: i
1471 REAL :: this_windspeed ( nlayer )
1472 REAL :: PI
1473 PARAMETER ( PI = 3.14159265359 )
1475 ! -------------------------------------------------------------------- !
1476 ! -------------------------------------------------------------------- !
1477 !~ Initialize variables
1478 ! --------------------
1479 dynamic = REAL ( 0 )
1481 !~ Compute the windspeed
1482 ! ---------------------
1483 DO i = 1, nlayer
1484 this_windspeed ( i ) = SQRT ( u(i)**2 + v(i)**2 )
1485 END DO
1487 !~ Compute the dynamic term
1488 ! -------------------------
1489 dynamic = ( this_windspeed(1)+this_windspeed(nlayer) ) / REAL (20)
1490 IF ( dynamic > REAL (2) ) dynamic = REAL ( 2 )
1491 dynamic = ( dynamic + REAL (3) ) / REAL ( 2 )
1492 dynamic = SIN ( dynamic*PI )
1493 dynamic = ( dynamic + REAL (1) ) / REAL ( 2 )
1496 END FUNCTION LLT_Windspeed
1499 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1500 !~ ~!
1501 !~ Name: ~!
1502 !~ LLT_Thermodynamic ~!
1503 !~ ~!
1504 !~ Description: ~!
1505 !~ This function computes the thermodynamic term for the low-level ~!
1506 !~ turbulence algorithm. ~!
1507 !~ ~!
1508 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1509 FUNCTION LLT_Thermodynamic ( nlayer, tK, hgt ) RESULT ( thermodynamic )
1510 IMPLICIT NONE
1512 !~ Variable Declaration
1513 ! --------------------
1514 INTEGER, INTENT ( IN ) :: nlayer
1515 REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperature (K)
1516 REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
1517 REAL :: thermodynamic
1519 !~ Internal function variables
1520 ! ---------------------------
1521 INTEGER :: i
1522 REAL :: lapse
1523 REAL :: PI
1524 PARAMETER ( PI = 3.14159265359 )
1526 ! -------------------------------------------------------------------- !
1527 ! -------------------------------------------------------------------- !
1529 !~ Initialize variables
1530 ! --------------------
1531 thermodynamic = REAL ( 0 )
1533 !~ Compute the lapse rate over the layer. The sign gets goofy here,
1534 !~ but works as coded below.
1535 ! -----------------------------------------------------------------
1536 lapse = ( tk(1) - tk(nlayer) ) * REAL ( 1000 )
1537 lapse = lapse / ( hgt(nlayer) - hgt(1) )
1539 !~ Compute the thermodynamic component
1540 ! -----------------------------------
1541 thermodynamic = lapse / REAL ( 10 )
1542 thermodynamic = ( thermodynamic + REAL (3) ) / REAL ( 2 )
1543 thermodynamic = SIN ( thermodynamic * PI )
1544 thermodynamic = ( thermodynamic + REAL (1) ) / REAL ( 2 )
1546 END FUNCTION LLT_Thermodynamic
1549 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1550 !~ ~!
1551 !~ Name: ~!
1552 !~ LLT_MountainWave ~!
1553 !~ ~!
1554 !~ Description: ~!
1555 !~ This function computes the mountain wave term for the low-level ~!
1556 !~ turbulence algorithm. ~!
1557 !~ ~!
1558 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1559 FUNCTION LLT_MountainWave ( nlayer, tdx, tdy, u, v, tK, hgt) &
1560 RESULT ( MountainWave )
1561 IMPLICIT NONE
1563 !~ Variable Declaration
1564 ! --------------------
1565 INTEGER, INTENT ( IN ) :: nlayer
1566 REAL, INTENT ( IN ) :: tdx !~ Terrain dx
1567 REAL, INTENT ( IN ) :: tdy !~ Terrain dy
1568 REAL, INTENT ( IN ) :: u ( nlayer ) !~ U components f
1569 REAL, INTENT ( IN ) :: v ( nlayer ) !~ V components
1570 REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperatures (K)
1571 REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
1572 REAL :: MountainWave !~ Mountain wave term
1574 !~ Internal function variables
1575 ! ---------------------------
1576 REAL :: u_term
1577 REAL :: v_term
1578 REAL :: uv_term
1579 REAL :: lapse
1580 REAL :: total_mw, this_total_mw
1581 REAL :: this_uv_term
1582 REAL :: min_uv_term, cross_terrain, max_total_mw
1583 INTEGER :: i, j, k
1585 REAL :: PI
1586 PARAMETER ( PI = 3.14159265359 )
1588 ! -------------------------------------------------------------------- !
1589 ! -------------------------------------------------------------------- !
1591 !~ Initialize variables
1592 ! --------------------
1593 MountainWave = REAL ( 0 )
1595 !~ Loop through the layer
1596 ! ----------------------
1597 DO i = 2, nlayer
1599 !~ Wind terrain term
1600 ! -----------------
1601 u_term = ( (u(i-1) + u(i) ) / REAL(2) ) * tdx
1602 v_term = ( (v(i-1) + v(i) ) / REAL(2) ) * tdy
1603 this_uv_term = ( u_term + v_term ) * REAL ( -1 )
1604 !IF ( uv_term < REAL (0) ) uv_term = REAL ( 0 )
1605 IF ( min_uv_term < REAL (0) ) min_uv_term = REAL ( 0 )
1606 IF ( i == 2 ) THEN
1607 !uv_term = this_uv_term
1608 min_uv_term = this_uv_term
1609 ELSE
1610 !IF ( this_uv_term < uv_term ) uv_term = this_uv_term
1611 IF ( this_uv_term < min_uv_term ) min_uv_term = this_uv_term
1612 END IF
1614 !~ Lapse rate
1615 ! ----------
1616 lapse = ( tK (i-1) - tK (i) ) * REAL ( 1000 )
1617 lapse = lapse / ABS ( hgt(i)-hgt(i-1) )
1618 IF ( lapse > REAL (0) ) lapse = REAL ( 0 )
1619 lapse = lapse * REAL ( -1 )
1621 this_total_mw = this_uv_term * lapse * REAL ( 40000 )
1622 IF ( i == 2 ) THEN
1623 total_mw = this_total_mw
1624 ELSE
1625 IF ( this_total_mw > total_mw ) total_mw = this_total_mw
1626 END IF
1628 END DO
1630 !min_uv_term = uv_term
1631 cross_terrain = min_uv_term * REAL ( 500 )
1633 IF ( min_uv_term < 0.03 ) THEN
1634 cross_terrain = REAL ( 0 )
1635 END IF
1637 IF ( cross_terrain > REAL (50) ) cross_terrain = REAL ( 50 )
1639 !~ Multiply the lapse (inversion) array and the mountain wave array
1640 ! ----------------------------------------------------------------
1641 IF ( total_mw > REAL (50) ) total_mw = REAL ( 50 )
1643 !~ Add the cross terrain flow and inversion term
1644 ! ---------------------------------------------
1645 MountainWave = ( total_mw*(cross_terrain/50.) ) + cross_terrain
1646 MountainWave = MountainWave / REAL ( 100 )
1648 END FUNCTION LLT_MountainWave
1652 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1653 !~ ~!
1654 !~ Name: ~!
1655 !~ LLT_TrappedWave ~!
1656 !~ ~!
1657 !~ Description: ~!
1658 !~ This function computes the trapped wave term for the low-level ~!
1659 !~ turbulence algorithm. ~!
1660 !~ ~!
1661 !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
1662 FUNCTION LLT_TrappedWave ( nlayer, u, v, p ) RESULT ( trapped )
1663 IMPLICIT NONE
1665 !~ Variable Declaration
1666 ! --------------------
1667 INTEGER, INTENT ( IN ) :: nlayer
1668 REAL, INTENT ( IN ) :: u ( nlayer )
1669 REAL, INTENT ( IN ) :: v ( nlayer )
1670 REAL, INTENT ( IN ) :: p ( nlayer )
1671 REAL :: trapped
1673 !~ Internal function variables
1674 ! ---------------------------
1675 INTEGER :: i
1676 REAL :: du, dv
1677 REAL :: scale_fact, this_p
1678 REAL :: dudv, this_dudv
1679 REAL :: PI
1680 PARAMETER ( PI = 3.14159265359 )
1682 !~ Scale parameters
1683 ! ----------------
1684 REAL, PARAMETER :: scale_950 = 0.050000 !~ 1/20
1685 REAL, PARAMETER :: scale_925 = 0.040000 !~ 1/25
1686 REAL, PARAMETER :: scale_900 = 0.025000 !~ 1/40
1687 REAL, PARAMETER :: scale_850 = 0.010000 !~ 1/100
1688 REAL, PARAMETER :: scale_800 = 0.005000 !~ 1/200
1689 REAL, PARAMETER :: scale_750 = 0.002941 !~ 1/340
1690 REAL, PARAMETER :: scale_700 = 0.001923 !~ 1/520
1691 REAL, PARAMETER :: scale_650 = 0.001351 !~ 1/740
1692 REAL, PARAMETER :: scale_600 = 0.001000 !~ 1/1000
1693 REAL, PARAMETER :: scale_550 = 0.000800 !~ 1/1250
1695 ! -------------------------------------------------------------------- !
1696 ! -------------------------------------------------------------------- !
1698 !~ Initialize variables
1699 ! --------------------
1700 trapped = REAL ( 0 )
1702 !~ Compute the trapped wave term
1703 ! ------------------
1704 dudv = REAL ( 0 )
1705 DO i = 2, nlayer
1707 !~ Compute dudv first
1708 ! ------------------
1709 du = u ( i-1 ) - u ( i )
1710 dv = v ( i-1 ) - v ( i )
1712 !~ Scale based on pressure level
1713 ! -----------------------------
1714 this_p = p ( i ) / REAL ( 100 )
1715 IF ( this_p > REAL (950) ) THEN
1716 scale_fact = scale_950
1717 ELSE IF ( this_p <= REAL (950) .AND. this_p > REAL (925) ) THEN
1718 scale_fact = scale_925
1719 ELSE IF ( this_p <= REAL (925) .AND. this_p > REAL (900) ) THEN
1720 scale_fact = scale_900
1721 ELSE IF ( this_p <= REAL (900) .AND. this_p > REAL (850) ) THEN
1722 scale_fact = scale_850
1723 ELSE IF ( this_p <= REAL (850) .AND. this_p > REAL (800) ) THEN
1724 scale_fact = scale_800
1725 ELSE IF ( this_p <= REAL (800) .AND. this_p > REAL (750) ) THEN
1726 scale_fact = scale_750
1727 ELSE IF ( this_p <= REAL (750) .AND. this_p > REAL (700) ) THEN
1728 scale_fact = scale_700
1729 ELSE IF ( this_p <= REAL (700) .AND. this_p > REAL (650) ) THEN
1730 scale_fact = scale_650
1731 ELSE IF ( this_p <= REAL (650) .AND. this_p > REAL (600) ) THEN
1732 scale_fact = scale_600
1733 ELSE IF ( this_p <= REAL (600) ) THEN
1734 scale_fact = scale_550
1735 END IF
1737 this_dudv = ( (du**2)*(dv**2) ) * scale_fact
1738 IF ( this_dudv > dudv ) dudv = this_dudv
1740 END DO
1742 trapped = dudv
1743 IF ( trapped > REAL ( 1 ) ) trapped = REAL ( 1 )
1744 trapped = trapped / REAL ( 4 )
1746 END FUNCTION LLT_TrappedWave
1748 END MODULE diag_functions
1749 #endif