#if (NMM_CORE == 1)
MODULE module_diag_afwa 1
CONTAINS
SUBROUTINE diag_afwa_stub
END SUBROUTINE diag_afwa_stub
END MODULE module_diag_afwa
#else
!WRF:MEDIATION_LAYER:PHYSICS
MODULE diag_functions 1
CONTAINS
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ calc_rh
!~
!~ Description:
!~ This function calculates relative humidity given pressure,
!~ temperature, and water vapor mixing ratio.
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_rh ( p, t, qv ) result ( rh ) 3
IMPLICIT NONE
REAL, INTENT(IN) :: p, t, qv
REAL :: rh
! Local
! -----
REAL, PARAMETER :: pq0=379.90516
REAL, PARAMETER :: a2=17.2693882
REAL, PARAMETER :: a3=273.16
REAL, PARAMETER :: a4=35.86
REAL, PARAMETER :: rhmin=1.
REAL :: q, qs
INTEGER :: i,j,k
! Following algorithms adapted from WRFPOST
! May want to substitute with another later
! -----------------------------------------
q=qv/(1.0+qv)
qs=pq0/p*exp(a2*(t-a3)/(t-a4))
rh=100.*q/qs
IF (rh .gt. 100.) THEN
rh=100.
ELSE IF (rh .lt. rhmin) THEN
rh=rhmin
ENDIF
END FUNCTION calc_rh
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ uv_wind
!~
!~ Description:
!~ This function calculates the wind speed given U and V wind
!~ components.
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION uv_wind ( u, v ) result ( wind_speed ) 3
IMPLICIT NONE
REAL, INTENT(IN) :: u, v
REAL :: wind_speed
wind_speed = sqrt( u*u + v*v )
END FUNCTION uv_wind
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ Theta
!~
!~ Description:
!~ This function calculates potential temperature as defined by
!~ Poisson's equation, given temperature and pressure ( hPa ).
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION Theta ( t, p ) 26
IMPLICIT NONE
!~ Variable declaration
! --------------------
REAL, INTENT ( IN ) :: t
REAL, INTENT ( IN ) :: p
REAL :: theta
REAL :: Rd ! Dry gas constant
REAL :: Cp ! Specific heat of dry air at constant pressure
REAL :: p0 ! Standard pressure ( 1000 hPa )
Rd = 287.04
Cp = 1004.67
p0 = 1000.00
!~ Poisson's equation
! ------------------
theta = t * ( (p0/p)**(Rd/Cp) )
END FUNCTION Theta
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ Thetae
!~
!~ Description:
!~ This function returns equivalent potential temperature using the
!~ method described in Bolton 1980, Monthly Weather Review, equation 43.
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION Thetae ( tK, p, rh, mixr ) 1
IMPLICIT NONE
!~ Variable Declarations
! ---------------------
REAL :: tK ! Temperature ( K )
REAL :: p ! Pressure ( hPa )
REAL :: rh ! Relative humidity
REAL :: mixr ! Mixing Ratio ( kg kg^-1)
REAL :: te ! Equivalent temperature ( K )
REAL :: thetae ! Equivalent potential temperature
REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
REAL, PARAMETER :: lv = 2.54*(10**6) ! Latent heat of vaporization
! (J kg^-1)
REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
! at pressure (J/deg kg)
REAL :: tlc ! LCL temperature
!~ Calculate the temperature of the LCL
! ------------------------------------
tlc = TLCL ( tK, rh )
!~ Calculate theta-e
! -----------------
thetae = (tK * (p0/p)**( (R/Cp)*(1.- ( (.28E-3)*mixr*1000.) ) ) )* &
exp( (((3.376/tlc)-.00254))*&
(mixr*1000.*(1.+(.81E-3)*mixr*1000.)) )
END FUNCTION Thetae
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ The2T.f90
!~
!~ Description:
!~ This function returns the temperature at any pressure level along a
!~ saturation adiabat by iteratively solving for it from the parcel
!~ thetae.
!~
!~ Dependencies:
!~ function thetae.f90
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
IMPLICIT NONE
!~ Variable Declaration
! --------------------
REAL, INTENT ( IN ) :: thetaeK
REAL, INTENT ( IN ) :: pres
LOGICAL, INTENT ( INOUT ) :: flag
REAL :: tparcel
REAL :: thetaK
REAL :: tovtheta
REAL :: tcheck
REAL :: svpr, svpr2
REAL :: smixr, smixr2
REAL :: thetae_check, thetae_check2
REAL :: tguess_2, correction
LOGICAL :: found
INTEGER :: iter
REAL :: R ! Dry gas constant
REAL :: Cp ! Specific heat for dry air
REAL :: kappa ! Rd / Cp
REAL :: Lv ! Latent heat of vaporization at 0 deg. C
R = 287.04
Cp = 1004.67
Kappa = R/Cp
Lv = 2.500E+6
!~ Make initial guess for temperature of the parcel
! ------------------------------------------------
tovtheta = (pres/100000.0)**(r/cp)
tparcel = thetaeK/exp(lv*.012/(cp*295.))*tovtheta
iter = 1
found = .false.
flag = .false.
DO
IF ( iter > 105 ) EXIT
tguess_2 = tparcel + REAL ( 1 )
svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
! ------------------------------------------------------------------ ~!
!~ When this function was orinially written, the final parcel ~!
!~ temperature check was based off of the parcel temperature and ~!
!~ not the theta-e it produced. As there are multiple temperature- ~!
!~ mixing ratio combinations that can produce a single theta-e value, ~!
!~ we change the check to be based off of the resultant theta-e ~!
!~ value. This seems to be the most accurate way of backing out ~!
!~ temperature from theta-e. ~!
!~ ~!
!~ Rentschler, April 2010 ~!
! ------------------------------------------------------------------ !
!~ Old way...
!thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
!tcheck = thetaK * tovtheta
!~ New way
thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
!~ Whew doggies - that there is some accuracy...
!IF ( ABS (tparcel-tcheck) < .05) THEN
IF ( ABS (thetaeK-thetae_check) < .001) THEN
found = .true.
flag = .true.
EXIT
END IF
!~ Old
!tparcel = tparcel + (tcheck - tparcel)*.3
!~ New
correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
tparcel = tparcel + correction
iter = iter + 1
END DO
!IF ( .not. found ) THEN
! print*, "Warning! Thetae to temperature calculation did not converge!"
! print*, "Thetae ", thetaeK, "Pressure ", pres
!END IF
END FUNCTION The2T
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ VirtualTemperature
!~
!~ Description:
!~ This function returns virtual temperature given temperature ( K )
!~ and mixing ratio.
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
IMPLICIT NONE
!~ Variable declaration
real, intent ( in ) :: tK !~ Temperature
real, intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
real :: Tv !~ Virtual temperature
Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
END FUNCTION VirtualTemperature
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ SaturationMixingRatio
!~
!~ Description:
!~ This function calculates saturation mixing ratio given the
!~ temperature ( K ) and the ambient pressure ( Pa ). Uses
!~ approximation of saturation vapor pressure.
!~
!~ References:
!~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
IMPLICIT NONE
REAL, INTENT ( IN ) :: tK
REAL, INTENT ( IN ) :: p
REAL :: ws
REAL :: es
es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
ws = ( 0.622*es ) / ( (p/100.0)-es )
END FUNCTION SaturationMixingRatio
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ tlcl
!~
!~ Description:
!~ This function calculates the temperature of a parcel of air would have
!~ if lifed dry adiabatically to it's lifting condensation level (lcl).
!~
!~ References:
!~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION TLCL ( tk, rh )
IMPLICIT NONE
REAL, INTENT ( IN ) :: tK !~ Temperature ( K )
REAL, INTENT ( IN ) :: rh !~ Relative Humidity ( % )
REAL :: tlcl
REAL :: denom, term1, term2
term1 = 1.0 / ( tK - 55.0 )
IF ( rh > REAL (0) ) THEN
term2 = ( LOG (rh/100.0) / 2840.0 )
ELSE
term2 = ( LOG (0.001/1.0) / 2840.0 )
END IF
denom = term1 - term2
tlcl = ( 1.0 / denom ) + REAL ( 55 )
END FUNCTION TLCL
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ PWat
!~
!~ Description:
!~ This function calculates precipitable water by summing up the
!~ water vapor in a column from the first eta layer to model top
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION Pwat ( nz, qv, qc, dz8w, rho ) 1
IMPLICIT NONE
!~ Variable declaration
! --------------------
INTEGER, INTENT ( IN ) :: nz !~ Number of vertical levels
REAL, INTENT ( IN ) :: qv ( nz ) !~ Specific humidity in layer (kg/kg)
REAL, INTENT ( IN ) :: qc ( nz ) !~ Cloud water in layer (kg/kg)
REAL, INTENT ( IN ) :: dz8w ( nz ) !~ Dist between vertical levels (m)
REAL, INTENT ( IN ) :: rho ( nz ) !~ Air density (kg/m^3)
REAL :: Pwat !~ Precipitable water (kg/m^2)
INTEGER :: k !~ Vertical index
!~ Precipitable water (kg/m^2)
! ---------------------------
Pwat=0
DO k = 1, nz
Pwat = Pwat + (qv(k) + qc(k)) * dz8w(k) * rho(k)
ENDDO
END FUNCTION Pwat
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ Buoyancy ~!
!~ ~!
!~ Description: ~!
!~ This function computes Convective Available Potential Energy (CAPE) ~!
!~ with inhibition as a result of water loading given the data required ~!
!~ to run up a sounding. ~!
!~ ~!
!~ Additionally, since we are running up a sounding anyways, this ~!
!~ function returns the height of the Level of Free Convection (LFC) and ~!
!~ the pressure at the LFC. That-a-ways, we don't have to run up a ~!
!~ sounding later, saving a relatively computationally expensive ~!
!~ routine. ~!
!~ ~!
!~ Usage: ~!
!~ ostat = Buoyancy ( tK, rh, p, hgt, sfc, CAPE, ZLFC, PLFC, parcel ) ~!
!~ ~!
!~ Where: ~!
!~ ~!
!~ IN ~!
!~ -- ~!
!~ tK = Temperature ( K ) ~!
!~ rh = Relative Humidity ( % ) ~!
!~ p = Pressure ( Pa ) ~!
!~ hgt = Geopotential heights ( m ) ~!
!~ sfc = integer rank within submitted arrays that represents the ~!
!~ surface ~!
!~ ~!
!~ OUT ~!
!~ --- ~!
!~ ostat INTEGER return status. Nonzero is bad. ~!
!~ CAPE ( J/kg ) Convective Available Potential Energy ~!
!~ ZLFC ( gpm ) Height at the LFC ~!
!~ PLFC ( Pa ) Pressure at the LFC ~!
!~ ~!
!~ tK, rh, p, and hgt are all REAL arrays, arranged from lower levels ~!
!~ to higher levels. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
parcel ) result (ostat)
IMPLICIT NONE
INTEGER, INTENT ( IN ) :: nz !~ Number of vertical levels
INTEGER, INTENT ( IN ) :: sfc !~ Surface level in the profile
REAL, INTENT ( IN ) :: tk ( nz ) !~ Temperature profile ( K )
REAL, INTENT ( IN ) :: rh ( nz ) !~ Relative Humidity profile ( % )
REAL, INTENT ( IN ) :: p ( nz ) !~ Pressure profile ( Pa )
REAL, INTENT ( IN ) :: hgt ( nz ) !~ Height profile ( gpm )
REAL, INTENT ( OUT ) :: cape !~ CAPE ( J kg^-1 )
REAL, INTENT ( OUT ) :: cin !~ CIN ( J kg^-1 )
REAL, INTENT ( OUT ) :: zlfc !~ LFC Height ( gpm )
REAL, INTENT ( OUT ) :: plfc !~ LFC Pressure ( Pa )
REAL, INTENT ( OUT ) :: lidx !~ Lifted index
INTEGER :: ostat !~ Function return status
!~ Nonzero is bad.
INTEGER, INTENT ( IN ) :: parcel !~ Most Unstable = 1 (default)
!~ Mean layer = 2
!~ Surface based = 3
!~ Derived profile variables
! -------------------------
REAL :: ws ( nz ) !~ Saturation mixing ratio
REAL :: w ( nz ) !~ Mixing ratio
REAL :: dTvK ( nz ) !~ Parcel / ambient Tv difference
REAL :: buoy ( nz ) !~ Buoyancy
REAL :: tlclK !~ LCL temperature ( K )
REAL :: plcl !~ LCL pressure ( Pa )
REAL :: nbuoy !~ Negative buoyancy
REAL :: pbuoy !~ Positive buoyancy
!~ Source parcel information
! -------------------------
REAL :: srctK !~ Source parcel temperature ( K )
REAL :: srcrh !~ Source parcel rh ( % )
REAL :: srcws !~ Source parcel sat. mixing ratio
REAL :: srcw !~ Source parcel mixing ratio
REAL :: srcp !~ Source parcel pressure ( Pa )
REAL :: srctheta !~ Source parcel theta ( K )
REAL :: srcthetaeK !~ Source parcel theta-e ( K )
INTEGER :: srclev !~ Level of the source parcel
REAL :: spdiff !~ Pressure difference
!~ Parcel variables
! ----------------
REAL :: ptK !~ Parcel temperature ( K )
REAL :: ptvK !~ Parcel virtual temperature ( K )
REAL :: tvK !~ Ambient virtual temperature ( K )
REAL :: pw !~ Parcel mixing ratio
!~ Other utility variables
! -----------------------
INTEGER :: i, j, k !~ Dummy iterator
INTEGER :: lfclev !~ Level of LFC
INTEGER :: prcl !~ Internal parcel type indicator
INTEGER :: mlev !~ Level for ML calculation
INTEGER :: lyrcnt !~ Number of layers in mean layer
LOGICAL :: flag !~ Dummy flag
LOGICAL :: wflag !~ Saturation flag
REAL :: freeze !~ Water loading multiplier
REAL :: pdiff !~ Pressure difference between levs
REAL :: pm, pu, pd !~ Middle, upper, lower pressures
REAL :: lidxu !~ Lifted index at upper level
REAL :: lidxd !~ Lifted index at lower level
!~ Thermo / dynamical constants
! ----------------------------
REAL :: Rd !~ Dry gas constant
PARAMETER ( Rd = 287.058 ) !~ J deg^-1 kg^-1
REAL :: Cp !~ Specific heat constant pressure
PARAMETER ( Cp = 1004.67 ) !~ J deg^-1 kg^-1
REAL :: g !~ Acceleration due to gravity
PARAMETER ( g = 9.80665 ) !~ m s^-2
REAL :: RUNDEF
PARAMETER ( RUNDEF = -9.999E30 )
!~ Initialize variables
! --------------------
ostat = 0
CAPE = REAL ( 0 )
CIN = REAL ( 0 )
ZLFC = RUNDEF
PLFC = RUNDEF
!~ Look for submitted parcel definition
!~ 1 = Most unstable
!~ 2 = Mean layer
!~ 3 = Surface based
! -------------------------------------
IF ( parcel > 3 .or. parcel < 1 ) THEN
prcl = 1
ELSE
prcl = parcel
END IF
!~ Initalize our parcel to be (sort of) surface based. Because of
!~ issues we've been observing in the WRF model, specifically with
!~ excessive surface moisture values at the surface, using a true
!~ surface based parcel is resulting a more unstable environment
!~ than is actually occuring. To address this, our surface parcel
!~ is now going to be defined as the parcel between 25-50 hPa
!~ above the surface. UPDATE - now that this routine is in WRF,
!~ going to trust surface info. GAC 20140415
! ----------------------------------------------------------------
!~ Compute mixing ratio values for the layer
! -----------------------------------------
DO k = sfc, nz
ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
w ( k ) = ( rh(k)/100.0 ) * ws ( k )
END DO
srclev = sfc
srctK = tK ( sfc )
srcrh = rh ( sfc )
srcp = p ( sfc )
srcws = ws ( sfc )
srcw = w ( sfc )
srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
!~ Compute the profile mixing ratio. If the parcel is the MU parcel,
!~ define our parcel to be the most unstable parcel within the lowest
!~ 180 mb.
! -------------------------------------------------------------------
mlev = sfc + 1
DO k = sfc + 1, nz
!~ Identify the last layer within 100 hPa of the surface
! -----------------------------------------------------
pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
IF ( pdiff <= REAL (100) ) mlev = k
!~ If we've made it past the lowest 180 hPa, exit the loop
! -------------------------------------------------------
IF ( pdiff >= REAL (180) ) EXIT
IF ( prcl == 1 ) THEN
!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
IF ( (w(k) > srcw) ) THEN
srctheta = Theta ( tK(k), p(k)/100.0 )
srcw = w ( k )
srclev = k
srctK = tK ( k )
srcrh = rh ( k )
srcp = p ( k )
END IF
END IF
END DO
!~ If we want the mean layer parcel, compute the mean values in the
!~ lowest 100 hPa.
! ----------------------------------------------------------------
lyrcnt = mlev - sfc + 1
IF ( prcl == 2 ) THEN
srclev = sfc
srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
srctheta = Theta ( srctK, srcp/100. )
END IF
srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
!~ Calculate temperature and pressure of the LCL
! ---------------------------------------------
tlclK = TLCL ( tK(srclev), rh(srclev) )
plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
!~ Now lift the parcel
! -------------------
buoy = REAL ( 0 )
pw = srcw
wflag = .false.
DO k = srclev, nz
IF ( p (k) <= plcl ) THEN
!~ The first level after we pass the LCL, we're still going to
!~ lift the parcel dry adiabatically, as we haven't added the
!~ the required code to switch between the dry adiabatic and moist
!~ adiabatic cooling. Since the dry version results in a greater
!~ temperature loss, doing that for the first step so we don't over
!~ guesstimate the instability.
! ----------------------------------------------------------------
IF ( wflag ) THEN
flag = .false.
!~ Above the LCL, our parcel is now undergoing moist adiabatic
!~ cooling. Because of the latent heating being undergone as
!~ the parcel rises above the LFC, must iterative solve for the
!~ parcel temperature using equivalant potential temperature,
!~ which is conserved during both dry adiabatic and
!~ pseudoadiabatic displacements.
! --------------------------------------------------------------
ptK = The2T ( srcthetaeK, p(k), flag )
!~ Calculate the parcel mixing ratio, which is now changing
!~ as we condense moisture out of the parcel, and is equivalent
!~ to the saturation mixing ratio, since we are, in theory, at
!~ saturation.
! ------------------------------------------------------------
pw = SaturationMixingRatio ( ptK, p(k) )
!~ Now we can calculate the virtual temperature of the parcel
!~ and the surrounding environment to assess the buoyancy.
! ----------------------------------------------------------
ptvK = VirtualTemperature ( ptK, pw )
tvK = VirtualTemperature ( tK (k), w (k) )
!~ Modification to account for water loading
! -----------------------------------------
freeze = 0.033 * ( 263.15 - pTvK )
IF ( freeze > 1.0 ) freeze = 1.0
IF ( freeze < 0.0 ) freeze = 0.0
!~ Approximate how much of the water vapor has condensed out
!~ of the parcel at this level
! ---------------------------------------------------------
freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
dTvK ( k ) = ptvK - tvK
buoy ( k ) = g * ( dTvK ( k ) / tvK )
ELSE
!~ Since the theta remains constant whilst undergoing dry
!~ adiabatic processes, can back out the parcel temperature
!~ from potential temperature below the LCL
! --------------------------------------------------------
ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
!~ Grab the parcel virtual temperture, can use the source
!~ mixing ratio since we are undergoing dry adiabatic cooling
! ----------------------------------------------------------
ptvK = VirtualTemperature ( ptK, srcw )
!~ Virtual temperature of the environment
! --------------------------------------
tvK = VirtualTemperature ( tK (k), w (k) )
!~ Buoyancy at this level
! ----------------------
dTvK ( k ) = ptvK - tvK
buoy ( k ) = g * ( dtvK ( k ) / tvK )
wflag = .true.
END IF
ELSE
!~ Since the theta remains constant whilst undergoing dry
!~ adiabatic processes, can back out the parcel temperature
!~ from potential temperature below the LCL
! --------------------------------------------------------
ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
!~ Grab the parcel virtual temperture, can use the source
!~ mixing ratio since we are undergoing dry adiabatic cooling
! ----------------------------------------------------------
ptvK = VirtualTemperature ( ptK, srcw )
!~ Virtual temperature of the environment
! --------------------------------------
tvK = VirtualTemperature ( tK (k), w (k) )
!~ Buoyancy at this level
! ---------------------
dTvK ( k ) = ptvK - tvK
buoy ( k ) = g * ( dtvK ( k ) / tvK )
END IF
!~ Chirp
! -----
! WRITE ( *,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw*1000.,ptvK,tvK,buoy(k)
END DO
!~ Add up the buoyancies, find the LFC
! -----------------------------------
flag = .false.
lfclev = -1
nbuoy = REAL ( 0 )
pbuoy = REAL ( 0 )
DO k = sfc + 1, nz
IF ( tK (k) < 253.15 ) EXIT
CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
!~ If we've already passed the LFC
! -------------------------------
IF ( flag .and. buoy (k) > REAL (0) ) THEN
pbuoy = pbuoy + buoy (k)
END IF
!~ We are buoyant now - passed the LFC
! -----------------------------------
IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
flag = .true.
pbuoy = pbuoy + buoy (k)
lfclev = k
END IF
!~ If we think we've passed the LFC, but encounter a negative layer
!~ start adding it up.
! ----------------------------------------------------------------
IF ( flag .and. buoy (k) < REAL (0) ) THEN
nbuoy = nbuoy + buoy (k)
!~ If the accumulated negative buoyancy is greater than the
!~ positive buoyancy, then we are capped off. Got to go higher
!~ to find the LFC. Reset positive and negative buoyancy summations
! ----------------------------------------------------------------
IF ( ABS (nbuoy) > pbuoy ) THEN
flag = .false.
nbuoy = REAL ( 0 )
pbuoy = REAL ( 0 )
lfclev = -1
END IF
END IF
END DO
!~ Calculate lifted index by interpolating difference between
!~ parcel and ambient Tv to 500mb.
! ----------------------------------------------------------
DO k = sfc + 1, nz
pm = 50000.
pu = p ( k )
pd = p ( k - 1 )
!~ If we're already above 500mb just set lifted index to 0.
!~ --------------------------------------------------------
IF ( pd .le. pm ) THEN
lidx = 0.
EXIT
ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
!~ Found trapping pressure: up, middle, down.
!~ We are doing first order interpolation.
! ------------------------------------------
lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
EXIT
ENDIF
END DO
!~ Assuming the the LFC is at a pressure level for now
! ---------------------------------------------------
IF ( lfclev > 0 ) THEN
PLFC = p ( lfclev )
ZLFC = hgt ( lfclev )
END IF
IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN
PLFC = REAL ( -1 )
ZLFC = REAL ( -1 )
END IF
IF ( CAPE /= CAPE ) cape = REAL ( 0 )
IF ( CIN /= CIN ) cin = REAL ( 0 )
!~ Chirp
! -----
! WRITE ( *,* ) ' CAPE: ', cape, ' CIN: ', cin
! WRITE ( *,* ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
! WRITE ( *,* ) ''
! WRITE ( *,* ) ' Exiting buoyancy.'
! WRITE ( *,* ) ' ==================================== '
! WRITE ( *,* ) ''
END FUNCTION Buoyancy
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: NGMSLP NMC SEA LEVEL PRESSURE REDUCTION
! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-02
!
! ABSTRACT:
!
! THIS ROUTINE COMPUTES SEA LEVEL PRESSURE USING THE
! HYDROSTATIC EQUATION WITH THE SHUELL CORRECTION. THE
! FOLLOWING IS BASED ON DOCUMENTATION IN SUBROUTINE
! OUTHYDRO OF THE NGM:
!
! THE FUNDAMENTAL HYDROSTATIC EQUATION IS
! D(HEIGHT)
! --------- = TAU = VIRTUAL TEMPERATURE * (RGAS/GRAVITY)
! D (Z)
! WHERE
! Z = MINUS LOG OF PRESSURE (-LN(P)).
!
! SEA-LEVEL PRESSURE IS COMPUTED FROM THE FORMULA
! PRESS(MSL) = PRESS(GROUND) * EXP( F)
! WHERE
! F = HEIGHT OF GROUND / MEAN TAU
! MEAN TAU = ( TAU(GRND) + TAU(SL) ) / 2
!
! IN THE NGM TAU(GRND) AND TAU(SL) ARE FIRST SET USING A
! 6.5DEG/KM LAPSE RATE FROM LOWEST MDL LEVEL. THIS IS MODIFIED
! BY A CORRECTION BASED ON THE CRITICAL TAU OF THE SHUELL
! CORRECTION:
! TAUCR=(RGASD/GRAVITY) * 290.66
!
! 1) WHERE ONLY TAU(SL) EXCEEDS TAUCR, CHANGE TAU(SL) TO TAUCR.
!
! 2) WHERE BOTH TAU(SL) AND TAU(GRND) EXCEED TAUCR,
! CHANGE TAU(SL) TO TAUCR-CONST*(TAU(GRND)-TAUCR )**2
! WHERE CONST = .005 (GRAVITY/RGASD)
!
! THE AVERAGE OF TAU(SL) AND TAU(GRND) IS THEN USED TOGETHER
! WITH THE GROUND HEIGHT AND PRESSURE TO DERIVE THE PRESSURE
! AT SEA LEVEL.
!
! HEIGHT OF THE 1000MB SURFACE IS COMPUTED FROM THE MSL PRESSURE
! FIELD USING THE FORMULA:
!
! P(MSL) - P(1000MB) = MEAN DENSITY * GRAVITY * HGT(1000MBS)
!
! WHERE P(MSL) IS THE SEA LEVEL PRESSURE FIELD WE HAVE JUST
! COMPUTED.
!
!
! MEB 6/13/02: THIS CODE HAS BEEN SIMPLIFIED CONSIDERABLY FROM
! THE ONE USED IN ETAPOST. HORIZONTAL SMOOTHING HAS BEEN
! REMOVED AND THE FIRST MODEL LEVEL IS USED RATHER
! THAN THE MEAN OF THE VIRTUAL TEMPERATURES IN
! THE LOWEST 30MB ABOVE GROUND TO COMPUTE TAU(GRND).
!
! .
!
! PROGRAM HISTORY LOG:
! 93-02-02 RUSS TREADON
! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D
! 00-01-04 JIM TUCCILLO - MPI VERSION
! 01-10-25 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
! 01-11-02 H CHUANG - MODIFIED LINE 234 FOR COMPUTATION OF
! SIGMA/HYBRID SLP
! 01-12-18 H CHUANG - INCLUDED SMOOTHING ALONG BOUNDARIES TO BE
! CONSISTENT WITH MESINGER SLP
! 02-06-13 MIKE BALDWIN - WRF VERSION
! 06-12-18 H CHUANG - BUG FIX TO CORRECT TAU AT SFC
! 14-04-17 G CREIGHTON - MODIFIED TO INSERT INTO AFWA DIAGNOSTICS IN WRF
!
!$$$
FUNCTION MSLP ( zsfc, psfc, zlev1, qlev1, tlev1 ) 6
implicit none
! DECLARE VARIABLES
REAL, INTENT ( IN ) :: zsfc !~ Surface height ( m )
REAL, INTENT ( IN ) :: psfc !~ Surface height ( m )
REAL, INTENT ( IN ) :: zlev1 !~ Level 1 height ( m )
REAL, INTENT ( IN ) :: qlev1 !~ Level 1 mixing ratio ( kg/kg )
REAL, INTENT ( IN ) :: tlev1 !~ Level 1 temperature ( K )
real,PARAMETER :: G=9.81
real,PARAMETER :: GI=1./G
real,PARAMETER :: RD=287.0
real,PARAMETER :: ZSL=0.0
real,PARAMETER :: TAUCR=RD*GI*290.66,CONST=0.005*G/RD
real,PARAMETER :: GORD=G/RD,DP=60.E2
real,PARAMETER :: GAMMA=6.5E-3
real MSLP,TVRT,TVRSFC,TAUSFC,TVRSL,TAUSL,TAUAVG
!
!**********************************************************************
! START NGMSLP HERE.
!
MSLP = PSFC
!
! COMPUTE LAYER TAU (VIRTUAL TEMP*RD/G).
TVRT = TLEV1*(1.0+0.608*QLEV1)
!TAU = TVRT*RD*GI
!
! COMPUTE TAU AT THE GROUND (Z=ZSFC) AND SEA LEVEL (Z=0)
! ASSUMING A CONSTANT LAPSE RATE OF GAMMA=6.5DEG/KM.
TVRSFC = TVRT + (ZLEV1 - ZSFC)*GAMMA
TAUSFC = TVRSFC*RD*GI
TVRSL = TVRT + (ZLEV1 - ZSL)*GAMMA
TAUSL = TVRSL*RD*GI
!
! IF NEED BE APPLY SHEULL CORRECTION.
IF ((TAUSL.GT.TAUCR).AND.(TAUSFC.LE.TAUCR)) THEN
TAUSL=TAUCR
ELSEIF ((TAUSL.GT.TAUCR).AND.(TAUSFC.GT.TAUCR)) THEN
TAUSL = TAUCR-CONST*(TAUSFC-TAUCR)**2
ENDIF
!
! COMPUTE MEAN TAU.
TAUAVG = 0.5*(TAUSL+TAUSFC)
!
! COMPUTE SEA LEVEL PRESSURE.
MSLP = PSFC*EXP(ZSFC/TAUAVG)
END FUNCTION MSLP
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ calc_fits ~!
!~ ~!
!~ Description: ~!
!~ This function computes Fighter Index Thermal Stress values given ~!
!~ dry bulb temperature, relative humidity, and pressure. ~!
!~ ~!
!~ Usage: ~!
!~ fitsval = calc_fits ( p, tK, rh ) ~!
!~ ~!
!~ Where: ~!
!~ p = Pressure ( Pa ) ~!
!~ tK = Temperature ( K ) ~!
!~ rh = Relative Humidity ( % ) ~!
!~ ~!
!~ Reference: ~!
!~ Stribley, R.F., S. Nunneley, 1978: Fighter Index of Thermal Stress: ~!
!~ Development of interim guidance for hot-weather USAF operations. ~!
!~ SAM-TR-78-6. Eqn. 9 ~!
!~ ~!
!~ Formula: ~!
!~ FITS = 0.8281*Twb + 0.3549*Tdb + 5.08 (degrees Celsius) ~!
!~ ~!
!~ Where: ~!
!~ Twb = Wet Bulb Temperature ~!
!~ Tdb = Dry Bulb Temperature ~!
!~ ~!
!~ Written: ~!
!~ Scott Rentschler, Software Engineering Services ~!
!~ Fine Scale Models Team ~!
!~ Air Force Weather Agency, 16WS/WXN ~!
!~ DSN: 271-3331 Comm: (402) 294-3331 ~!
!~ scott.rentschler@offutt.af.mil ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_fits ( p, tK, rh ) RESULT ( fits )
implicit none
!~ Variable declaration
! --------------------
real, intent ( in ) :: p !~ Pressure ( Pa )
real, intent ( in ) :: tK !~ Temperature ( K )
real, intent ( in ) :: rh !~ Rel Humidity ( % )
real :: fits !~ FITS index value
!~ Utility variables
! --------------------------
real :: twb !~ Wet bulb temperature ( C )
real :: wbt
! ---------------------------------------------------------------------- !
! ---------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
fits = REAL ( 0 )
!~ Get the wet bulb temperature in degrees Celsius
! -----------------------------------------------
twb = WetBulbTemp ( p, tK, rh ) - 273.15
!~ Compute the FITS index
! ----------------------
fits = 0.8281*twb + 0.3549*( tK - 273.15 ) + 5.08
!~ Convert the index to Kelvin
! ---------------------------
fits = fits + 273.15
END FUNCTION calc_fits
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ calc_wc
!~
!~ Description:
!~ This function calculates wind chill given temperature ( K ) and
!~ wind speed ( m/s )
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_wc ( tK, wspd ) RESULT ( wc )
implicit none
!~ Variable Declarations
! ---------------------
real, intent ( in ) :: tK
real, intent ( in ) :: wspd
real :: tF, wc, wspd_mph
wspd_mph = wspd * 2.23693629 ! convert to mph
tF = (( tK - 273.15 ) * ( REAL (9) / REAL (5) ) ) + REAL ( 32 )
wc = 35.74 &
+ ( 0.6215 * tF ) &
- ( 35.75 * ( wspd_mph**0.16 ) ) &
+ ( 0.4275 * tF * ( wspd_mph**0.16 ) )
wc = (( wc - REAL (32) ) * ( REAL (5) / REAL (9) ) ) + 273.15
END FUNCTION calc_wc
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ calc_hi
!~
!~ Description:
!~ This subroutine calculates the heat index. Requires temperature ( K )
!~ and relative humidity ( % ).
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_hi ( Tk, RH ) result ( HI )
implicit none
!~ Variable declarations
! ---------------------
real, intent ( in ) :: Tk
real, intent ( in ) :: RH
real :: tF, tF2, rh2, HI
!~ If temperature > 70F then calculate heat index, else set it equal
!~ to dry temperature
! -----------------------------------------------------------------
IF ( Tk > 294.26111 ) THEN
tF = ( (Tk - 273.15) * (REAL (9)/REAL (5)) ) + REAL ( 32 )
tF2 = tF ** 2
rh2 = RH ** 2
HI = -42.379 &
+ ( 2.04901523 * tF ) &
+ ( 10.14333127 * RH ) &
- ( 0.22475541 * tF * RH ) &
- ( 6.83783E-03 * tF2 ) &
- ( 5.481717E-02 * rh2 ) &
+ ( 1.22874E-03 * tF2 * RH ) &
+ ( 8.5282E-04 * tF * rh2 ) &
- ( 1.99E-06 * tF2 * rh2 )
HI = ((HI - REAL (32)) * (REAL (5)/REAL (9))) + 273.15
ELSE
HI = Tk
END IF
END FUNCTION calc_hi
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ WetBulbTemp ~!
!~ ~!
!~ Description: ~!
!~ This function approximates the Wet Bulb Temperature (K) provided ~!
!~ dry bulb temperature (K), relative humidity (%), and pressure (Pa). ~!
!~ ~!
!~ Usage: ~!
!~ wbt = WetBulbTemperature ( p, tK, rh ) ~!
!~ ~!
!~ Where: ~!
!~ p = Pressure ( Pa ) ~!
!~ tK = Temperature ( K ) ~!
!~ rh = Relative Humidity ( % ) ~!
!~ ~!
!~ Reference: ~!
!~ American Society of Civil Engineers ~!
!~ Evapotraspiration and Irrigation Water Requirements ~!
!~ Jensen et al (1990) ASCE Manual No. 70, pp 176-177 ~!
!~ ~!
!~ Written: ~!
!~ Scott Rentschler, Software Engineering Services ~!
!~ Fine Scale Models Team ~!
!~ Air Force Weather Agency ~!
!~ DSM: 271-3331 Comm: (402) 294-3331 ~!
!~ scott.rentschler@offutt.af.mil ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION WetBulbTemp ( p, tK, rh) result( wbt )
implicit none
!~ Variable delclaration
! ---------------------
real, intent ( in ) :: p !~ Pressure ( Pa )
real, intent ( in ) :: tK !~ Temperature ( K )
real, intent ( in ) :: rh !~ Relative Humidity ( % )
real :: wbt !~ Wet Bulb Temperature ( K )
!~ Utility variables
! -----------------
real :: tdK !~ Dewpoint temperature ( K )
real :: tC !~ Temperature ( C )
real :: tdC !~ Dewpoint temperature ( K )
real :: svapr !~ Saturation vapor pressure ( Pa )
real :: vapr !~ Ambient vapor pressure ( Pa )
real :: gamma !~ Dummy term
real :: delta !~ Dummy term
! ---------------------------------------------------------------------- !
! ---------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
wbt = REAL ( 0 )
tC = tK - 273.15
!~ Compute saturation vapor pressure ( Pa )
! ----------------------------------------
svapr = calc_es ( tK ) * REAL ( 100 )
!~ Compute vapor pressure
! ----------------------
vapr = svapr * ( rh / REAL (100) )
!~ Grab the dewpoint
! -----------------
tdC = calc_Dewpoint ( tC, rh )
tdK = tdC + 273.15
!~ Compute dummy terms
! -------------------
gamma = 0.00066 * ( p / REAL (1000) )
delta = REAL ( 4098 ) * ( vapr / REAL(1000) ) / ( (tC+237.3)**2 )
!~ Compute the wet bulb temperature
! --------------------------------
wbt = ( ((gamma * tC) + (delta * tdC)) / (gamma + delta) ) + 273.15
END FUNCTION WetBulbTemp
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ calc_Dewpoint
!~
!~ Description:
!~ This function approximates dewpoint given temperature and rh.
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_Dewpoint ( tC, rh) result( Dewpoint )
implicit none
!~ Variable Declaration
! --------------------
real, intent ( in ) :: tC
real, intent ( in ) :: rh
real :: Dewpoint
real :: term, es, e1, e, logs, expon
expon = ( 7.5*tC ) / ( 237.7+tC )
es = 6.112 * ( 10**expon ) ! Saturated vapor pressure
e = es * ( rh/100.0 ) ! Vapor pressure
logs = LOG10 ( e/6.112 )
Dewpoint = ( 237.7*logs ) / ( 7.5-logs )
END FUNCTION calc_Dewpoint
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~
!~ Name:
!~ calc_es
!~
!~ Description:
!~ This function returns the saturation vapor pressure over water ( hPa )
!~ given temperature ( K ).
!~
!~ References:
!~ The algorithm is due to Nordquist, W.S., 1973: "Numerical approximations
!~ of selected meteorological parameters for cloud physics problems,"
!~ ecom-5475, Atmospheric Sciences Laboratory, U.S. Army Electronics
!~ Command, White Sands Missile Range, New Mexico, 88002
!~
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION calc_es ( tK ) result ( es )
implicit none
!~ Variable Declaration
! --------------------
real, intent ( in ) :: tK
real :: es
real :: p1, p2, c1
p1 = 11.344 - 0.0303998 * tK
p2 = 3.49149 - 1302.8844 / tK
c1 = 23.832241 - 5.02808 * ALOG10 ( tK )
es = 10.**(c1-1.3816E-7*10.**p1+8.1328E-3*10.**p2-2949.076/tK)
END FUNCTION calc_es
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ CATTurbulence ~!
!~ ~!
!~ Description: ~!
!~ This function calculates the turbulence index ( TI ) for one layer ~!
!~ in the atmosphere given the horizontal wind components and the geo- ~!
!~ potential height of two pressure levels. The index is computed for ~!
!~ the layer between the levels using the deformation and convergence ~!
!~ of the wind field at the top and bottom of the layer and the vertical ~!
!~ wind shear is calculated within the layer. The equation used for ~!
!~ calculating TI is given by: ~!
!~ ~!
!~ ~!
!~ TI = VWS * ( DEF + CONV ) ~!
!~ ~!
!~ Where: ~!
!~ VWS = Vertical wind shear ~!
!~ DEF = Deformation ~!
!~ CONV = Convergence ~!
!~ ~!
!~ Notes: ~!
!~ ~!
!~ References: ~!
!~ Ellrod, G.P. and D.J. Knapp, An objective clear-air turbulence ~!
!~ forecasting technique: verification and operational use, Weather ~!
!~ and Forecasting, 7, March 1992, pp. 150-165. ~!
!~ ~!
!~ Written: ~!
!~ Scott Rentschler, Software Engineering Services ~!
!~ Fine Scale Models Team ~!
!~ Air Force Weather Agency, 16WS/WXN ~!
!~ DSN: 271-3331 Comm: (402) 294-3331 ~!
!~ scott.rentschler@offutt.af.mil ~!
!~ ~!
!~ History: ~!
!~ 1 February 2008 ................... Scott Rentschler, (SES), 2WG/WEA ~!
!~ INITIAL VERSION ~!
!~ ~!
!~ 8 July 2009 ....................... Scott Rentschler, (SES), 2WG/WEA ~!
!~ Adapted for new driver. ~!
!~ ~!
!~ 1 November 2012 ......................... Scott Rentschler, 16WS/WXN ~!
!~ Modified to accept layer argument, which adds the flexibility to make ~!
!~ the computation for whichever flight level is desired. Cleaned up ~!
!~ some of the code and added a couple comments. ~!
!~ ~!
!~ 28 August 2013 .................... Braedi Wickard, SEMS/NG/16WS/WXN ~!
!~ Adapted for use within the Unified Post Processor. UPP can not handle ~!
!~ the layer argument for flight levels, so reverted to hardcoded levels ~!
!~ ~!
!~ 25 April 2014 ............................. Glenn Creighton, 16WS/WXN ~!
!~ Adapted for use within WRF. WRF already computes many of these terms. ~!
!~ Stripped everything down to its bare bones to remove need to compute ~!
!~ horizontal terms, now using deformation variables already within WRF. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION CATTurbulence ( ugrdbot, ugrdtop, vgrdbot, vgrdtop &
,defor11bot, defor11top, defor12bot, defor12top &
,defor22bot, defor22top, zbot, ztop ) result ( ti )
IMPLICIT NONE
!~ Variable declarations
! ---------------------
REAL, INTENT ( IN ) :: ugrdbot !~ U-wind bottom of layer
REAL, INTENT ( IN ) :: ugrdtop !~ U-wind top of layer
REAL, INTENT ( IN ) :: vgrdbot !~ V-wind bottom of layer
REAL, INTENT ( IN ) :: vgrdtop !~ V-wind top of layer
REAL, INTENT ( IN ) :: defor11bot !~ 2*du/dx at bottom of layer
REAL, INTENT ( IN ) :: defor11top !~ 2*du/dx at top of layer
REAL, INTENT ( IN ) :: defor12bot !~ du/dy+dv/dx at bottom of layer
REAL, INTENT ( IN ) :: defor12top !~ du/dy+dv/dx at top of layer
REAL, INTENT ( IN ) :: defor22bot !~ 2*dv/dy at bottom of layer
REAL, INTENT ( IN ) :: defor22top !~ 2*dv/dy at top of layer
REAL, INTENT ( IN ) :: zbot !~ Height grid bottom
REAL, INTENT ( IN ) :: ztop !~ Height grid top
REAL :: ti !~ Turbulence index
!~ Function utility variables
! --------------------------
REAL :: dudx, dudx1, dudx2 !~ Wind differentials
REAL :: dvdy, dvdy1, dvdy2
REAL :: dudz, dvdz
REAL :: depth, vws, conv !~ Depth, vertical wind shear, convergence
REAL :: def, shear, stretch !~ Deformation, shear, stretching terms
!~ Initialize variables.
! ----------------------
ti = REAL ( 0 )
!~ Compute vertical wind shear
! ---------------------------
depth = ABS ( ztop - zbot )
dudz = ( ugrdbot - ugrdtop ) / depth
dvdz = ( vgrdbot - vgrdtop ) / depth
vws = SQRT ( dudz**2 + dvdz**2 )
dudx1 = defor11top / 2.
dudx2 = defor11bot / 2.
dudx = ( dudx1 + dudx2 ) / REAL ( 2 )
dvdy1 = defor22top / 2.
dvdy2 = defor22bot / 2.
dvdy = ( dvdy1 + dvdy2 ) / REAL ( 2 )
!~ Compute the deformation
! -----------------------
stretch = dudx - dvdy
shear = ( defor12top + defor12bot ) / REAL ( 2 )
def = SQRT ( stretch**2 + shear**2 )
!~ Compute the convergence
! -----------------------
conv = - ( dudx + dvdy )
!~ Compute the turbulence index
! ----------------------------
ti = vws * ( def + conv ) * 1.0E+07
IF ( ti /= ti ) ti = REAL ( 0 )
IF ( ti < 0 ) ti = REAL ( 0 )
END FUNCTION CATTurbulence
FUNCTION lin_interp ( x, f, y ) result ( g ),4
! Purpose:
! interpolates f(x) to point y
! assuming f(x) = f
(x0) + a * (x - x0)
! where a = ( f(x1) - f(x0) ) / (x1 - x0)
! x0 <= x <= x1
! assumes x is monotonically increasing
! Author: D. Fillmore :: J. Done changed from r8 to r4
! Pilfered for AFWA diagnostics - G Creighton
implicit none
real, intent(in), dimension(:) :: x ! grid points
real, intent(in), dimension(:) :: f ! grid function values
real, intent(in) :: y ! interpolation point
real :: g ! interpolated function value
integer :: k ! interpolation point index
integer :: n ! length of x
real :: a
n = size(x)
! find k such that x(k) < y =< x(k+1)
! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
if (y <= x(1)) then
k = 1
else if (y >= x(n)) then
k = n - 1
else
k = 1
do while (y > x(k+1) .and. k < n)
k = k + 1
end do
end if
! interpolate
a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
g = f(k) + a * (y - x(k))
END FUNCTION lin_interp
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ LLT_Windspeed ~!
!~ ~!
!~ Description: ~!
!~ This function computes the dynamic term for the low-level turbulence ~!
!~ algorithm. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION LLT_Windspeed ( nlayer, u, v ) RESULT ( dynamic )
IMPLICIT NONE
!~ Variable Declaration
! --------------------
INTEGER, INTENT ( IN ) :: nlayer
REAL, INTENT ( IN ) :: u ( nlayer )
REAL, INTENT ( IN ) :: v ( nlayer )
REAL :: dynamic
!~ Internal function variables
! ---------------------------
INTEGER :: i
REAL :: this_windspeed ( nlayer )
REAL :: PI
PARAMETER ( PI = 3.14159265359 )
! -------------------------------------------------------------------- !
! -------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
dynamic = REAL ( 0 )
!~ Compute the windspeed
! ---------------------
DO i = 1, nlayer
this_windspeed ( i ) = SQRT ( u(i)**2 + v(i)**2 )
END DO
!~ Compute the dynamic term
! -------------------------
dynamic = ( this_windspeed(1)+this_windspeed(nlayer) ) / REAL (20)
IF ( dynamic > REAL (2) ) dynamic = REAL ( 2 )
dynamic = ( dynamic + REAL (3) ) / REAL ( 2 )
dynamic = SIN ( dynamic*PI )
dynamic = ( dynamic + REAL (1) ) / REAL ( 2 )
END FUNCTION LLT_Windspeed
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ LLT_Thermodynamic ~!
!~ ~!
!~ Description: ~!
!~ This function computes the thermodynamic term for the low-level ~!
!~ turbulence algorithm. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION LLT_Thermodynamic ( nlayer, tK, hgt ) RESULT ( thermodynamic )
IMPLICIT NONE
!~ Variable Declaration
! --------------------
INTEGER, INTENT ( IN ) :: nlayer
REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperature (K)
REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
REAL :: thermodynamic
!~ Internal function variables
! ---------------------------
INTEGER :: i
REAL :: lapse
REAL :: PI
PARAMETER ( PI = 3.14159265359 )
! -------------------------------------------------------------------- !
! -------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
thermodynamic = REAL ( 0 )
!~ Compute the lapse rate over the layer. The sign gets goofy here,
!~ but works as coded below.
! -----------------------------------------------------------------
lapse = ( tk(1) - tk(nlayer) ) * REAL ( 1000 )
lapse = lapse / ( hgt(nlayer) - hgt(1) )
!~ Compute the thermodynamic component
! -----------------------------------
thermodynamic = lapse / REAL ( 10 )
thermodynamic = ( thermodynamic + REAL (3) ) / REAL ( 2 )
thermodynamic = SIN ( thermodynamic * PI )
thermodynamic = ( thermodynamic + REAL (1) ) / REAL ( 2 )
END FUNCTION LLT_Thermodynamic
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ LLT_MountainWave ~!
!~ ~!
!~ Description: ~!
!~ This function computes the mountain wave term for the low-level ~!
!~ turbulence algorithm. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION LLT_MountainWave ( nlayer, tdx, tdy, u, v, tK, hgt) &
RESULT ( MountainWave )
IMPLICIT NONE
!~ Variable Declaration
! --------------------
INTEGER, INTENT ( IN ) :: nlayer
REAL, INTENT ( IN ) :: tdx !~ Terrain dx
REAL, INTENT ( IN ) :: tdy !~ Terrain dy
REAL, INTENT ( IN ) :: u ( nlayer ) !~ U components f
REAL, INTENT ( IN ) :: v ( nlayer ) !~ V components
REAL, INTENT ( IN ) :: tK ( nlayer ) !~ Temperatures (K)
REAL, INTENT ( IN ) :: hgt ( nlayer ) !~ Heights ( m )
REAL :: MountainWave !~ Mountain wave term
!~ Internal function variables
! ---------------------------
REAL :: u_term
REAL :: v_term
REAL :: uv_term
REAL :: lapse
REAL :: total_mw, this_total_mw
REAL :: this_uv_term
REAL :: min_uv_term, cross_terrain, max_total_mw
INTEGER :: i, j, k
REAL :: PI
PARAMETER ( PI = 3.14159265359 )
! -------------------------------------------------------------------- !
! -------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
MountainWave = REAL ( 0 )
!~ Loop through the layer
! ----------------------
DO i = 2, nlayer
!~ Wind terrain term
! -----------------
u_term = ( (u(i-1) + u(i) ) / REAL(2) ) * tdx
v_term = ( (v(i-1) + v(i) ) / REAL(2) ) * tdy
this_uv_term = ( u_term + v_term ) * REAL ( -1 )
!IF ( uv_term < REAL (0) ) uv_term = REAL ( 0 )
IF ( min_uv_term < REAL (0) ) min_uv_term = REAL ( 0 )
IF ( i == 2 ) THEN
!uv_term = this_uv_term
min_uv_term = this_uv_term
ELSE
!IF ( this_uv_term < uv_term ) uv_term = this_uv_term
IF ( this_uv_term < min_uv_term ) min_uv_term = this_uv_term
END IF
!~ Lapse rate
! ----------
lapse = ( tK (i-1) - tK (i) ) * REAL ( 1000 )
lapse = lapse / ABS ( hgt(i)-hgt(i-1) )
IF ( lapse > REAL (0) ) lapse = REAL ( 0 )
lapse = lapse * REAL ( -1 )
this_total_mw = this_uv_term * lapse * REAL ( 40000 )
IF ( i == 2 ) THEN
total_mw = this_total_mw
ELSE
IF ( this_total_mw > total_mw ) total_mw = this_total_mw
END IF
END DO
!min_uv_term = uv_term
cross_terrain = min_uv_term * REAL ( 500 )
IF ( min_uv_term < 0.03 ) THEN
cross_terrain = REAL ( 0 )
END IF
IF ( cross_terrain > REAL (50) ) cross_terrain = REAL ( 50 )
!~ Multiply the lapse (inversion) array and the mountain wave array
! ----------------------------------------------------------------
IF ( total_mw > REAL (50) ) total_mw = REAL ( 50 )
!~ Add the cross terrain flow and inversion term
! ---------------------------------------------
MountainWave = ( total_mw*(cross_terrain/50.) ) + cross_terrain
MountainWave = MountainWave / REAL ( 100 )
END FUNCTION LLT_MountainWave
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
!~ ~!
!~ Name: ~!
!~ LLT_TrappedWave ~!
!~ ~!
!~ Description: ~!
!~ This function computes the trapped wave term for the low-level ~!
!~ turbulence algorithm. ~!
!~ ~!
!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
FUNCTION LLT_TrappedWave ( nlayer, u, v, p ) RESULT ( trapped ),1
IMPLICIT NONE
!~ Variable Declaration
! --------------------
INTEGER, INTENT ( IN ) :: nlayer
REAL, INTENT ( IN ) :: u ( nlayer )
REAL, INTENT ( IN ) :: v ( nlayer )
REAL, INTENT ( IN ) :: p ( nlayer )
REAL :: trapped
!~ Internal function variables
! ---------------------------
INTEGER :: i
REAL :: du, dv
REAL :: scale_fact, this_p
REAL :: dudv, this_dudv
REAL :: PI
PARAMETER ( PI = 3.14159265359 )
!~ Scale parameters
! ----------------
REAL, PARAMETER :: scale_950 = 0.050000 !~ 1/20
REAL, PARAMETER :: scale_925 = 0.040000 !~ 1/25
REAL, PARAMETER :: scale_900 = 0.025000 !~ 1/40
REAL, PARAMETER :: scale_850 = 0.010000 !~ 1/100
REAL, PARAMETER :: scale_800 = 0.005000 !~ 1/200
REAL, PARAMETER :: scale_750 = 0.002941 !~ 1/340
REAL, PARAMETER :: scale_700 = 0.001923 !~ 1/520
REAL, PARAMETER :: scale_650 = 0.001351 !~ 1/740
REAL, PARAMETER :: scale_600 = 0.001000 !~ 1/1000
REAL, PARAMETER :: scale_550 = 0.000800 !~ 1/1250
! -------------------------------------------------------------------- !
! -------------------------------------------------------------------- !
!~ Initialize variables
! --------------------
trapped = REAL ( 0 )
!~ Compute the trapped wave term
! ------------------
dudv = REAL ( 0 )
DO i = 2, nlayer
!~ Compute dudv first
! ------------------
du = u ( i-1 ) - u ( i )
dv = v ( i-1 ) - v ( i )
!~ Scale based on pressure level
! -----------------------------
this_p = p ( i ) / REAL ( 100 )
IF ( this_p > REAL (950) ) THEN
scale_fact = scale_950
ELSE IF ( this_p <= REAL (950) .AND. this_p > REAL (925) ) THEN
scale_fact = scale_925
ELSE IF ( this_p <= REAL (925) .AND. this_p > REAL (900) ) THEN
scale_fact = scale_900
ELSE IF ( this_p <= REAL (900) .AND. this_p > REAL (850) ) THEN
scale_fact = scale_850
ELSE IF ( this_p <= REAL (850) .AND. this_p > REAL (800) ) THEN
scale_fact = scale_800
ELSE IF ( this_p <= REAL (800) .AND. this_p > REAL (750) ) THEN
scale_fact = scale_750
ELSE IF ( this_p <= REAL (750) .AND. this_p > REAL (700) ) THEN
scale_fact = scale_700
ELSE IF ( this_p <= REAL (700) .AND. this_p > REAL (650) ) THEN
scale_fact = scale_650
ELSE IF ( this_p <= REAL (650) .AND. this_p > REAL (600) ) THEN
scale_fact = scale_600
ELSE IF ( this_p <= REAL (600) ) THEN
scale_fact = scale_550
END IF
this_dudv = ( (du**2)*(dv**2) ) * scale_fact
IF ( this_dudv > dudv ) dudv = this_dudv
END DO
trapped = dudv
IF ( trapped > REAL ( 1 ) ) trapped = REAL ( 1 )
trapped = trapped / REAL ( 4 )
END FUNCTION LLT_TrappedWave
END MODULE diag_functions
MODULE module_diag_afwa 1
USE diag_functions
CONTAINS
SUBROUTINE afwa_diagnostics_driver ( grid , config_flags & 1,34
, moist &
, scalar &
, chem &
, th_phy , pi_phy , p_phy &
, u_phy , v_phy &
, dz8w , p8w , t8w , rho_phy , rho &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe &
, its, ite, jts, jte &
, k_start, k_end )
USE module_domain, ONLY : domain , domain_clock_get
USE module_configure, ONLY : grid_config_rec_type, model_config_rec
USE module_state_description
USE module_model_constants
USE module_utility
USE module_streams, ONLY: history_alarm, auxhist2_alarm
#ifdef DM_PARALLEL
USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval
#endif
IMPLICIT NONE
TYPE ( domain ), INTENT(INOUT) :: grid
TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
INTEGER :: k_start , k_end, its, ite, jts, jte
REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_moist), &
INTENT(IN ) :: moist
REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_scalar), &
INTENT(IN ) :: scalar
REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_chem), &
INTENT(IN ) :: chem
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: th_phy &
, pi_phy &
, p_phy &
, u_phy &
, v_phy &
, dz8w &
, p8w &
, t8w &
, rho_phy &
, rho
! Local
! -----
CHARACTER*512 :: message
CHARACTER*256 :: timestr
INTEGER :: i,j,k,nz,ostat
INTEGER :: icing_opt
REAL :: bdump
INTEGER :: i_start, i_end, j_start, j_end
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qrain &
, qsnow &
, qgrpl &
, qvapr &
, qcloud &
, qice &
, ncloud &
, ngraup &
, rh &
, rh_cld &
, ptot &
, z_e &
, zagl
REAL, DIMENSION( ims:ime, jms:jme, 5 ) :: dustc
REAL, DIMENSION( ims:ime, jms:jme ) :: rh2m &
, rh20m &
, tv2m &
, tv20m &
, wind10m &
, wup_mask &
, wind125m &
, llws &
, pwater
LOGICAL :: do_buoy_calc
REAL :: zlfc_msl, dum1, dum2, dum3, wind_vel, wind_blend
REAL :: prate_mm_per_hr, factor
REAL :: u1km, v1km, ublend, vblend, u2000, v2000, us, vs
LOGICAL :: is_target_level
! Timing
! ------
TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime
TYPE(WRFU_TimeInterval) :: dtint, histint, aux2int
LOGICAL :: is_after_history_dump, is_output_timestep, is_first_timestep
! Chirp the routine name for debugging purposes
! ---------------------------------------------
write ( message, * ) 'inside afwa_diagnostics_driver'
CALL wrf_debug( 100 , message )
! Get timing info
! Want to know if when the last history output was
! Check history and auxhist2 alarms to check last ring time and how often
! they are set to ring
! -----------------------------------------------------------------------
CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, &
ringinterval=histint)
CALL WRFU_ALARMGET( grid%alarms( AUXHIST2_ALARM ), prevringtime=aux2_time, &
ringinterval=aux2int)
! Get domain clock
! ----------------
CALL domain_clock_get ( grid, current_time=CurrTime, &
simulationStartTime=StartTime, &
current_timestr=timestr, time_step=dtint )
! Set some booleans for use later
! Following uses an overloaded .lt.
! ---------------------------------
is_after_history_dump = ( Currtime .lt. hist_time + dtint )
! Following uses an overloaded .ge.
! ---------------------------------
is_output_timestep = (Currtime .ge. hist_time + histint - dtint .or. &
Currtime .ge. aux2_time + aux2int - dtint )
write ( message, * ) 'is output timestep? ', is_output_timestep
CALL wrf_debug( 100 , message )
! Following uses an overloaded .eq.
! ---------------------------------
is_first_timestep = ( Currtime .eq. StartTime + dtint )
! Here is an optional check for bad data in case the model has gone
! off the hinges unchecked until now. This happens under certain
! configurations and can lead to very wrong data writes. This is just
! a simple check for sane winds and potential temperatures.
! --------------------------------------------------------------------
IF ( config_flags%afwa_bad_data_check .GT. 0 ) THEN
IF ( ( is_output_timestep ) .AND. ( .NOT. is_first_timestep ) ) THEN
DO i=its, MIN( ide-1, ite )
DO k=k_start, k_end
DO j=jts, MIN( jde-1, jte )
IF ( ( u_phy(i,k,j) .GT. 300. ) .OR. &
( u_phy(i,k,j) .LT. -300. ) .OR. &
( v_phy(i,k,j) .GT. 300. ) .OR. &
( v_phy(i,k,j) .LT. -300. ) .OR. &
( th_phy(i,k,j) .GT. 9999. ) .OR. &
( th_phy(i,k,j) .LT. 99. ) ) THEN
write ( message, * ) "AFWA Diagnostics: ERROR - Model winds and/or " // &
"potential temperature appear to be bad. If you do not want this check, " // &
"set afwa_bad_data_check=0. i=",i,", j=",j,", k=",k,", u_phy=",u_phy(i,k,j), &
", v_phy=", v_phy(i,k,j),", th_phy=",th_phy(i,k,j)
CALL wrf_error_fatal( message )
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
ENDIF
! 3-D arrays for moisture variables
! ---------------------------------
DO i=ims, ime
DO k=kms, kme
DO j=jms, jme
qvapr(i,k,j) = moist(i,k,j,P_QV)
qrain(i,k,j) = moist(i,k,j,P_QR)
qsnow(i,k,j) = moist(i,k,j,P_QS)
qgrpl(i,k,j) = moist(i,k,j,P_QG)
qcloud(i,k,j) = moist(i,k,j,P_QC)
qice(i,k,j) = moist(i,k,j,P_QI)
ncloud(i,k,j) = scalar(i,k,j,P_QNC)
ENDDO
ENDDO
ENDDO
! Total pressure
! --------------
DO i=ims, ime
DO k=kms, kme
DO j=jms, jme
ptot(i,k,j)=grid%pb(i,k,j)+grid%p(i,k,j)
ENDDO
ENDDO
ENDDO
! ZAGL (height above ground)
! --------------------------
DO i=ims, ime
DO k=kms, kme
DO j=jms, jme
zagl(i,k,j)=grid%z(i,k,j)-grid%ht(i,j)
ENDDO
ENDDO
ENDDO
! Calculate relative humidity
! ---------------------------
DO i=ims,ime
DO k=kms,kme
DO j=jms,jme
rh(i,k,j)=calc_rh(ptot(i,k,j),grid%t_phy(i,k,j), qvapr(i,k,j))
rh_cld(i,k,j)=calc_rh(ptot(i,k,j),grid%t_phy(i,k,j), qvapr(i,k,j)+qcloud(i,k,j)+qice(i,k,j))
ENDDO
ENDDO
ENDDO
! Time-step precipitation (convective + nonconvective)
! --------------------------------------------------------------
DO i=ims,ime
DO j=jms,jme
grid % afwa_precip(i,j) = grid%raincv(i,j) + grid%rainncv(i,j)
grid % afwa_totprecip(i,j) = grid%rainc(i,j) + grid%rainnc(i,j)
ENDDO
ENDDO
! Calculate precipitable water
! ----------------------------
nz=kme-kms+1
DO i=ims,ime
DO j=jms,jme
grid % afwa_pwat ( i, j ) = Pwat( nz, &
qvapr(i,kms:kme,j), &
qcloud(i,kms:kme,j), &
dz8w(i,kms:kme,j), &
rho(i,kms:kme,j) )
ENDDO
ENDDO
! After each history dump, reset max/min value arrays
! ----------------------------------------------------------------------
IF ( is_after_history_dump ) THEN
DO j = jms, jme
DO i = ims, ime
grid % wspd10max(i,j) = 0.
grid % afwa_llws(i,j) = 0.
ENDDO
ENDDO
ENDIF
! Calculate the max 10 m wind speed between output times
! ------------------------------------------------------
! UPDATE 20150112 - GAC
! Diagnose from model 10 m winds, and blend with 1 km AGL
! winds when precipitation rate is > 50 mm/hr to account
! for increased surface wind gust potential when precip
! is heavy and when winds aloft are strong. Will use the
! higher of the surface and the blended winds. Blending
! is linear weighted between 50-150 mm/hr precip rates.
! -------------------------------------------------------
DO j = jms, jme
DO i = ims, ime
wind_vel = uv_wind ( grid % u10(i,j) , grid % v10(i,j) )
prate_mm_per_hr = ( grid % afwa_precip(i,j) / grid % dt ) * 3600.
! Is this an area of heavy precip? Calculate 1km winds to blend down
! -------------------------------------------------------------------
IF ( prate_mm_per_hr .GT. 50. ) THEN
is_target_level=.false.
DO k=kms,kme
IF ( ( zagl(i,k,j) >= 1000. ) .and. &
( .NOT. is_target_level ) .and. &
( k .ne. kms ) ) THEN
is_target_level = .true.
u1km = u_phy(i,k-1,j) + (1000. - (zagl(i,k-1,j))) &
* ((u_phy(i,k,j) - u_phy(i,k-1,j))/(zagl(i,k,j)))
v1km = v_phy(i,k-1,j) + (1000. - (zagl(i,k-1,j))) &
* ((v_phy(i,k,j) - v_phy(i,k-1,j))/(zagl(i,k,j)))
EXIT ! We've found our level, break the loop
ENDIF
ENDDO
! Compute blended wind
! --------------------
factor = MAX ( ( ( 150. - prate_mm_per_hr ) / 100. ), 0. )
ublend = grid % u10(i,j) * factor + u1km * (1. - factor)
vblend = grid % v10(i,j) * factor + v1km * (1. - factor)
wind_blend = uv_wind ( ublend, vblend )
! Set the surface wind to the blended wind if higher
! --------------------------------------------------
IF ( wind_blend .GT. wind_vel ) THEN
wind_vel = wind_blend
ENDIF
ENDIF
IF ( wind_vel .GT. grid % wspd10max(i,j) ) THEN
grid % wspd10max(i,j) = wind_vel
ENDIF
ENDDO
ENDDO
! Calculate 0-2000 foot (0 - 609.6 meter) shear.
! ----------------------------------------------
DO j = jts, jte
DO i = its, ite
is_target_level=.false.
DO k=kms,kme
IF ( ( zagl(i,k,j) >= 609.6 ) .and. &
( .NOT. is_target_level ) .and. &
( k .ne. kms ) ) THEN
is_target_level = .true.
u2000 = u_phy(i,k-1,j) + (609.6 - (zagl(i,k-1,j))) &
* ((u_phy(i,k,j) - u_phy(i,k-1,j))/(zagl(i,k,j)))
v2000 = v_phy(i,k-1,j) + (609.6 - (zagl(i,k-1,j))) &
* ((v_phy(i,k,j) - v_phy(i,k-1,j))/(zagl(i,k,j)))
us = u2000 - grid % u10(i,j)
vs = v2000 - grid % v10(i,j)
llws(i,j) = uv_wind ( us , vs )
IF ( llws(i,j) .gt. grid % afwa_llws(i,j) ) THEN
grid % afwa_llws(i,j) = llws(i,j)
ENDIF
EXIT ! We've found our level, break the loop
ENDIF
ENDDO
ENDDO
ENDDO
#if ( WRF_CHEM == 1 )
! Surface dust concentration array (ug m-3)
! -----------------------------------------
DO i=ims, ime
DO j=jms, jme
dustc(i,j,1)=chem(i,k_start,j,p_dust_1)*rho(i,k_start,j)
dustc(i,j,2)=chem(i,k_start,j,p_dust_2)*rho(i,k_start,j)
dustc(i,j,3)=chem(i,k_start,j,p_dust_3)*rho(i,k_start,j)
dustc(i,j,4)=chem(i,k_start,j,p_dust_4)*rho(i,k_start,j)
dustc(i,j,5)=chem(i,k_start,j,p_dust_5)*rho(i,k_start,j)
ENDDO
ENDDO
#else
dustc(ims:ime,jms:jme,:)=0.
#endif
! Calculate severe weather diagnostics. These variables should only be
! output at highest frequency output. (e.g. auxhist2)
! ---------------------------------------------------------------------
IF ( config_flags % afwa_severe_opt == 1 ) THEN
! After each history dump, reset max/min value arrays
! Note: This resets up_heli_max which is currently calculated within
! rk_first_rk_step_part2.F, may want to move to this diagnostics package
! later
! ----------------------------------------------------------------------
IF ( is_after_history_dump ) THEN
DO j = jms, jme
DO i = ims, ime
! grid%wspd10max(i,j) = 0.
grid%w_up_max(i,j) = 0.
grid%w_dn_max(i,j) = 0.
grid%tcoli_max(i,j) = 0.
grid%grpl_flx_max(i,j) = 0.
grid%up_heli_max(i,j) = 0.
! grid%refd_max(i,j) = 0.
grid%afwa_tornado(i,j) = 0.
grid%midrh_min_old(i,j) = grid%midrh_min(i,j) ! Save old midrh_min
grid%midrh_min(i,j) = 999.
grid%afwa_hail(i,j) = 0.
ENDDO
ENDDO
ENDIF ! is_after_history_dump
IF ( ( is_first_timestep ) .OR. ( is_output_timestep ) ) THEN
do_buoy_calc = .true.
ELSE
do_buoy_calc = .false.
ENDIF
!-->RAS
! We need to do some neighboring gridpoint comparisons in this next function;
! set these values so we don't go off the edges of the domain. Updraft
! duration on domain edges will always be 0.
! ----------------------------------------------------------------------
i_start = its
i_end = ite
j_start = jts
j_end = jte
IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
config_flags%nested) i_start = MAX( ids+1, its )
IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
config_flags%nested) i_end = MIN( ide-1, ite )
IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
config_flags%nested) j_start = MAX( jds+1, jts )
IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
config_flags%nested) j_end = MIN( jde-1, jte )
IF ( config_flags%periodic_x ) i_start = its
IF ( config_flags%periodic_x ) i_end = ite
CALL severe_wx_diagnostics ( grid % wspd10max &
, grid % w_up_max &
, grid % w_dn_max &
, grid % up_heli_max &
, grid % tcoli_max &
, grid % midrh_min_old &
, grid % midrh_min &
, grid % afwa_hail &
, grid % afwa_cape &
, grid % afwa_cin &
! , grid % afwa_cape_mu &
! , grid % afwa_cin_mu &
, grid % afwa_zlfc &
, grid % afwa_plfc &
, grid % afwa_lidx &
, llws &
, grid % afwa_tornado &
, grid % grpl_flx_max &
, grid % u10 &
, grid % v10 &
, grid % w_2 &
, grid % uh &
, grid % t_phy &
, grid % t2 &
, grid % z &
, grid % ht &
, grid % tornado_mask &
, grid % tornado_dur &
, grid % dt &
, grid % afwa_pwat &
, u_phy &
, v_phy &
, ptot &
, qice &
, qsnow &
, qgrpl &
, ngraup &
, qvapr, qrain, qcloud &
, rho &
, dz8w &
, rh &
, do_buoy_calc &
, ims, ime, jms, jme, kms, kme &
, its, ite, jts, jte &
, k_start, k_end &
, j_start, j_end, i_start, i_end )
ENDIF ! afwa_severe_opt == 1
! Calculate precipitation type diagnostics
! ----------------------------------------
IF ( config_flags % afwa_ptype_opt == 1 ) THEN
! First initialize precip buckets
! -------------------------------
IF ( grid % itimestep .eq. 1) THEN
DO i=ims,ime
DO j=jms,jme
grid % afwa_rain(i,j)=0.
grid % afwa_snow(i,j)=0.
grid % afwa_ice(i,j)=0.
grid % afwa_fzra(i,j)=0.
grid % afwa_snowfall(i,j)=0.
ENDDO
ENDDO
ENDIF
! Diagnose precipitation type
! ---------------------------
CALL precip_type_diagnostics ( grid % t_phy &
, grid % t2 &
, rh &
, grid % z &
, grid % ht &
, grid % afwa_precip &
, grid % swdown &
, grid % afwa_rain &
, grid % afwa_snow &
, grid % afwa_ice &
, grid % afwa_fzra &
, grid % afwa_snowfall &
, grid % afwa_ptype_ccn_tmp &
, grid % afwa_ptype_tot_melt &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_ptype_opt == 1
! --------------------------------------------------------------
! The following packages are calculated only on output timesteps
! --------------------------------------------------------------
IF ( is_output_timestep ) THEN
! Calculate mean sea level pressure
! ---------------------------------
DO j = jms, jme
DO i = ims, ime
grid % afwa_mslp ( i, j ) = MSLP ( grid % ht ( i, j ) &
, grid % psfc ( i, j ) &
, grid % z ( i, kms, j ) &
, qvapr ( i, kms, j ) &
, grid % t_phy ( i, kms, j ) )
ENDDO
ENDDO
! Calculate 10 meter winds
! ------------------------
DO i=ims,ime
DO j=jms,jme
wind10m(i,j)=uv_wind(grid%u10(i,j),grid%v10(i,j))
ENDDO
ENDDO
! Calculate 2 meter relative humidity/Tv
! --------------------------------------
DO i=ims,ime
DO j=jms,jme
rh2m(i,j)=calc_rh(grid%psfc(i,j), grid%t2(i,j), grid%q2(i,j))
tv2m(i,j)=grid%t2(i,j) * (1 + 0.61 * grid%q2(i,j))
ENDDO
ENDDO
! Calculate buoyancy parameters.
! ------------------------------
IF ( config_flags % afwa_buoy_opt == 1 ) THEN
nz = k_end - k_start + 1
! Calculate buoyancy for surface parcel
! -------------------------------------
DO j = jts, jte
DO i = its, ite
ostat = Buoyancy ( nz &
, grid%t_phy(i,kms:kme,j) &
, rh(i,kms:kme,j) &
, ptot(i,kms:kme,j) &
, grid % z(i,kms:kme,j) &
, 1 &
, grid % afwa_cape(i,j) &
, grid % afwa_cin(i,j) &
, zlfc_msl &
, grid % afwa_plfc(i,j) &
, grid % afwa_lidx(i,j) &
, 3 ) !Surface
! Subtract terrain height to convert ZLFC from MSL to AGL
! -------------------------------------------------------
IF ( zlfc_msl .ge. grid % ht ( i, j ) ) THEN
grid % afwa_zlfc ( i, j ) = zlfc_msl - grid % ht ( i, j )
ELSE
grid % afwa_zlfc( i, j ) = -1.
ENDIF
! Add 273.15 to lifted index per some standard I don't understand
! ---------------------------------------------------------------
IF ( grid % afwa_lidx ( i, j ) .ne. 999. ) THEN
grid % afwa_lidx ( i, j ) = grid % afwa_lidx ( i, j ) + 273.15
ENDIF
! Calculate buoyancy again for most unstable layer
! ------------------------------------------------
ostat = Buoyancy ( nz &
, grid%t_phy(i,kms:kme,j) &
, rh(i,kms:kme,j) &
, ptot(i,kms:kme,j) &
, grid % z(i,kms:kme,j) &
, 1 &
, grid % afwa_cape_mu(i,j) &
, grid % afwa_cin_mu(i,j) &
, dum1 &
, dum2 &
, dum3 &
, 1 ) !Most unstable
ENDDO
ENDDO
ENDIF ! afwa_buoy_opt == 1
IF ( config_flags % afwa_therm_opt == 1 ) THEN
write ( message, * ) 'Calculating thermal indices'
CALL wrf_debug( 100 , message )
CALL thermal_diagnostics ( grid % t2 &
, grid % psfc &
, rh2m &
, wind10m &
, grid % afwa_heatidx &
, grid % afwa_wchill &
, grid % afwa_fits &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_therm_opt == 1
IF ( config_flags % afwa_turb_opt == 1 ) THEN
write ( message, * ) 'Calculating turbulence indices'
!~ For now, hard code turbulence layer bottom and top AGL heights
! --------------------------------------------------------------
grid % afwa_tlyrbot = (/ 1500., 3000., 4600., 6100., 7600., 9100., &
10700. /)
grid % afwa_tlyrtop = (/ 3000., 4600., 6100., 7600., 9100., 10700., &
12700. /)
call turbulence_diagnostics ( u_phy &
, v_phy &
, grid % t_phy &
, ptot &
, zagl &
, grid % defor11 &
, grid % defor12 &
, grid % defor22 &
, grid % afwa_turb &
, grid % afwa_llturb &
, grid % afwa_llturblgt &
, grid % afwa_llturbmdt &
, grid % afwa_llturbsvr &
!, config_flags % num_turb_layers &
, 7 &
, grid % afwa_tlyrbot &
, grid % afwa_tlyrtop &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_turb_opt == 1
! Calculate equivalent radar reflectivity factor (z_e) using
! old RIP code (2004) if running radar or VIL packages.
! ----------------------------------------------------------
IF ( config_flags % afwa_radar_opt == 1 .or. &
config_flags % afwa_vil_opt == 1 ) THEN
write ( message, * ) 'Calculating Radar'
CALL wrf_debug( 100 , message )
CALL wrf_dbzcalc ( rho &
, grid%t_phy &
, qrain &
, qsnow &
, qgrpl &
, z_e &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_radar_opt == 1 .or. afwa_vil_opt == 1
! Calculate derived radar variables
! ---------------------------------
! UPDATE: removed refd_max calculation, which was not working correctly.
! To correctly calculate refd_max, this routine could be called every
! time step, but instead, we are only going to calculate reflectivity
! on output time steps and avoid cost of calculating refd_max. GAC2014
! ----------------------------------------------------------------------
IF ( config_flags % afwa_radar_opt == 1 ) THEN
write ( message, * ) 'Calculating derived radar variables'
CALL wrf_debug( 100 , message )
CALL radar_diagnostics ( grid % refd &
, grid % refd_com &
! , grid % refd_max &
, grid % echotop &
, grid % z &
, z_e &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_radar_opt == 1
! Calculate VIL and reflectivity every history output timestep
! ------------------------------------------------------------
IF ( config_flags % afwa_vil_opt == 1 ) THEN
write ( message, * ) 'Calculating VIL'
CALL wrf_debug( 100 , message )
CALL vert_int_liquid_diagnostics ( grid % vil &
, grid % radarvil &
, grid % t_phy &
, qrain &
, qsnow &
, qgrpl &
, z_e &
, dz8w &
, rho &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_vil_opt ==1
! Calculate icing and freezing level
! ----------------------------------
IF ( config_flags % afwa_icing_opt == 1 ) THEN
! Determine icing option from microphysics scheme
! -----------------------------------------------
IF ( config_flags % mp_physics == GSFCGCESCHEME ) THEN
icing_opt=1
ELSEIF ( config_flags % mp_physics == ETAMPNEW ) THEN
icing_opt=2
ELSEIF ( config_flags % mp_physics == THOMPSON ) THEN
icing_opt=3
ELSEIF ( config_flags % mp_physics == WSM5SCHEME .OR. &
config_flags % mp_physics == WSM6SCHEME ) THEN
icing_opt=4
ELSEIF ( config_flags % mp_physics == MORR_TWO_MOMENT ) THEN
! Is this run with prognostic cloud droplets or no?
! -------------------------------------------------
IF (config_flags % progn > 0) THEN
icing_opt=6
ELSE
icing_opt=5
ENDIF
ELSEIF ( config_flags % mp_physics == WDM6SCHEME ) THEN
icing_opt=7
ELSE
icing_opt=0 ! Not supported
ENDIF
write ( message, * ) 'Calculating Icing with icing opt ',icing_opt
CALL wrf_debug( 100 , message )
CALL icing_diagnostics ( icing_opt &
, grid % fzlev &
, grid % icing_lg &
, grid % icing_sm &
, grid % qicing_lg_max &
, grid % qicing_sm_max &
, grid % qicing_lg &
, grid % qicing_sm &
, grid % icingtop &
, grid % icingbot &
, grid % t_phy &
, grid % z &
, dz8w &
, rho &
, qrain &
, qcloud &
, ncloud &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF ! afwa_icing_opt
! Calculate visiblility diagnostics
! ---------------------------------
IF ( config_flags % afwa_vis_opt == 1 ) THEN
! Interpolate 20m temperature and RH
! ----------------------------------
DO i=ims,ime
DO j=jms,jme
tv20m(i,j) = -999.
rh20m(i,j) = -999.
DO k = kps, MIN(kpe,kde-1)
IF (tv20m (i,j) .eq. -999. .AND. zagl (i,k,j) .ge. 20.) THEN
! If the lowest model level > 20 m AGL, interpolate using 2-m temperature/RH
! --------------------------------------------------------------------------
IF (k .eq. kps) THEN
tv20m(i,j) = tv2m(i,j) + &
(20. - 2.) / &
(zagl(i,k,j) - 2.) * &
(grid%t_phy(i,k,j) * (1 + 0.61 * qvapr(i,k,j)) - tv2m(i,j))
rh20m(i,j) = rh2m(i,j) + &
(20. - 2.) / &
(zagl(i,k,j) - 2.) * &
(rh(i,k,j) - rh2m(i,j))
ELSE
tv20m(i,j) = grid%t_phy(i,k-1,j) * (1 + 0.61 * qvapr(i,k-1,j)) + &
((20. - zagl(i,k-1,j)) / &
(zagl(i,k,j) - zagl(i,k-1,j))) * &
(grid%t_phy(i,k,j) * (1 + 0.61 * qvapr(i,k,j)) - &
grid%t_phy(i,k-1,j) * (1 + 0.61 * qvapr(i,k-1,j)))
rh20m(i,j) = rh (i,k-1,j) + &
((20. - zagl (i,k-1,j)) / &
(zagl (i,k,j) - zagl (i,k-1,j))) * &
(rh (i,k,j) - rh (i,k-1,j))
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
! Calculate 125 meter winds
! -------------------------
DO i=ims,ime
DO j=jms,jme
wind125m(i,j) = -999.
DO k = kps, MIN(kpe,kde-1)
IF (wind125m (i,j) .eq. -999. .AND. zagl (i,k,j) .ge. 125.) THEN
! If the lowest model level > 125 m AGL, use level 1 wind
! -------------------------------------------------------
IF (k .eq. kps) THEN
wind125m(i,j) = uv_wind(u_phy(i,k,j),v_phy(i,k,j))
ELSE
wind125m(i,j) = uv_wind(u_phy(i,k-1,j),v_phy(i,k-1,j)) + &
((125. - zagl(i,k-1,j)) / &
(zagl(i,k,j) - zagl(i,k-1,j))) * &
(uv_wind(u_phy(i,k,j),v_phy(i,k,j)) - &
uv_wind(u_phy(i,k-1,j),v_phy(i,k-1,j)))
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
write ( message, * ) 'Calculating visibility'
CALL wrf_debug( 100 , message )
CALL vis_diagnostics ( qcloud(ims:ime,k_start,jms:jme) &
, qrain(ims:ime,k_start,jms:jme) &
, qice(ims:ime,k_start,jms:jme) &
, qsnow(ims:ime,k_start,jms:jme) &
, qgrpl(ims:ime,k_start,jms:jme) &
, rho(ims:ime,k_start,jms:jme) &
, wind10m &
, wind125m &
, grid % afwa_pwat &
, grid % q2 &
, rh2m &
, rh20m &
, tv2m &
, tv20m &
, dustc &
, grid % afwa_vis &
, grid % afwa_vis_dust &
, grid % afwa_vis_alpha &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF
! Calculate cloud diagnostics
! ---------------------------
IF ( config_flags % afwa_cloud_opt == 1 ) THEN
write ( message, * ) 'Calculating cloud'
CALL wrf_debug( 100 , message )
CALL cloud_diagnostics (qcloud &
, qice &
, qsnow &
, rh_cld &
, dz8w &
, rho &
, grid % z &
, grid % ht &
, grid % afwa_cloud &
, grid % afwa_cloud_ceil &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
ENDIF
ENDIF ! is_output_timestep
END SUBROUTINE afwa_diagnostics_driver
SUBROUTINE severe_wx_diagnostics ( wspd10max & 1
, w_up_max &
, w_dn_max &
, up_heli_max &
, tcoli_max &
, midrh_min_old &
, midrh_min &
, afwa_hail &
, cape &
, cin &
! , cape_mu &
! , cin_mu &
, zlfc &
, plfc &
, lidx &
, llws &
, afwa_tornado &
, grpl_flx_max &
, u10 &
, v10 &
, w_2 &
, uh &
, t_phy &
, t2 &
, z &
, ht &
, tornado_mask &
, tornado_dur &
, dt &
, pwat &
, u_phy &
, v_phy &
, p &
, qi &
, qs &
, qg &
, ng &
, qv, qr, qc &
, rho &
, dz8w &
, rh &
, do_buoy_calc &
, ims, ime, jms, jme, kms, kme &
, its, ite, jts, jte &
, k_start, k_end &
, j_start, j_end, i_start, i_end )
INTEGER, INTENT(IN) :: its, ite, jts, jte, k_start, k_end &
, ims, ime, jms, jme, kms, kme &
, j_start, j_end, i_start, i_end
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: p &
, w_2 &
, t_phy &
, u_phy &
, v_phy &
, qi &
, qs &
, qg &
, ng &
, qv, qr, qc &
, rho &
, z &
, dz8w &
, rh
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: u10 &
, v10 &
, wspd10max &
, uh &
, t2 &
, ht &
, midrh_min_old &
, up_heli_max &
, llws &
, pwat
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: w_up_max &
, w_dn_max &
, tcoli_max &
, midrh_min &
, afwa_hail &
, afwa_tornado &
, grpl_flx_max &
, tornado_mask &
, tornado_dur
! REAL, DIMENSION( ims:ime, jms:jme ), &
! INTENT( OUT) :: cape &
! , cin &
! , cape_mu &
! , cin_mu &
! , zlfc &
! , plfc &
! , lidx
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: cape &
, cin &
, zlfc &
, plfc &
, lidx
REAL, INTENT(IN) :: dt
LOGICAL, INTENT(IN) :: do_buoy_calc
! Local
! -----
INTEGER :: i,j,k
INTEGER :: kts,kte
REAL :: zagl, zlfc_msl, melt_term, midrh_term, hail, midrh
REAL :: dum1, dum2, dum3
REAL :: tornado, lfc_term, shr_term, midrh2_term, uh_term
REAL :: wind_vel, p_tot, tcoli, grpl_flx, w_n15, qg_n15
INTEGER :: nz, ostat
REAL, DIMENSION( ims:ime, jms:jme ) :: w_up &
, w_dn &
, tornado_mask_prev &
, tornado_dur_prev
REAL :: time_factor, time_factor_prev
LOGICAL :: is_target_level
! Calculate midlevel relative humidity minimum
! --------------------------------------------
DO i=ims,ime
DO j=jms,jme
is_target_level=.false.
DO k=kms,kme
zagl = z(i,k,j) - ht(i,j)
IF ( ( zagl >= 3500. ) .and. &
( .NOT. is_target_level ) .and. &
( k .ne. kms ) ) THEN
is_target_level = .true.
midrh = rh(i,k-1,j) + (3500. - (z(i,k-1,j) - ht(i,j))) &
* ((rh(i,k,j) - rh(i,k-1,j))/(z(i,k,j) - z(i,k-1,j)))
IF ( midrh .lt. midrh_min(i,j) ) THEN
midrh_min(i,j) = midrh
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
! Calculate the max 10 m wind speed between output times
! ------------------------------------------------------
!DO j = jts, jte
! DO i = its, ite
! wind_vel = uv_wind ( u10(i,j) , v10(i,j) )
! IF ( wind_vel .GT. wspd10max(i,j) ) THEN
! wspd10max(i,j) = wind_vel
! ENDIF
! ENDDO
!ENDDO
! Vertical velocity quantities between output times
! -------------------------------------------------
w_up=0.
w_dn=0.
DO j = jts, jte
DO k = k_start, k_end
DO i = its, ite
p_tot = p(i,k,j) / 100.
! Check vertical velocity field below 400 mb
IF ( p_tot .GT. 400. .AND. w_2(i,k,j) .GT. w_up(i,j) ) THEN
w_up(i,j) = w_2(i,k,j)
IF ( w_up(i,j) .GT. w_up_max(i,j) ) THEN
w_up_max(i,j) = w_up(i,j)
ENDIF
ENDIF
IF ( p_tot .GT. 400. .AND. w_2(i,k,j) .LT. w_dn(i,j) ) THEN
w_dn(i,j) = w_2(i,k,j)
IF ( w_dn(i,j) .LT. w_dn_max(i,j) ) THEN
w_dn_max(i,j) = w_dn(i,j)
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
! Hail diameter in millimeters (Weibull distribution)
! ---------------------------------------------------
DO j = jts, jte
DO i = its, ite
melt_term=max(t2(i,j)-288.15,0.)
midrh_term=max(2*(min(midrh_min(i,j),midrh_min_old(i,j))-70.),0.)
! Change exponent to 1.1 to reduce probabilities for large sizes
!hail=max((w_up(i,j)/1.4)**1.25-melt_term-midrh_term,0.)
hail=max((w_up(i,j)/1.4)**1.1-melt_term-midrh_term,0.)
hail=hail*((uh(i,j)/100)+0.25)
IF ( hail .gt. afwa_hail(i,j) ) THEN
afwa_hail(i,j)=hail
ENDIF
ENDDO
ENDDO
! Lightning (total column-integrated cloud ice)
! Note this formula is basically stolen from the VIL calculation.
! ---------------------------------------------------------------
DO j = jts, jte
DO i = its, ite
tcoli=0.
DO k = k_start, k_end
tcoli = tcoli + &
(qi (i,k,j) + &
qs (i,k,j) + &
qg (i,k,j)) &
*dz8w (i,k,j) * rho(i,k,j)
ENDDO
IF ( tcoli .GT. tcoli_max(i,j) ) THEN
tcoli_max(i,j) = tcoli
ENDIF
ENDDO
ENDDO
! Lighting (pseudo graupel flux calculation)
! Model graupel mixing ration (g/kg) times w (m/s) at the -15C level
! Values should range from around 25 in marginal lightning situations,
! to over 400 in situations with very frequent lightning.
! --------------------------------------------------------------------
DO j = jts, jte
DO i = its, ite
grpl_flx=0
w_n15=-999.
DO k = k_start, k_end
! Interpolate w and qg to -15C level and calculate graupel flux
! as simply graupel x vertical velocity at -15C
! -------------------------------------------------------------
IF ( t_phy (i,k,j) .LE. 258.15 .AND. w_n15 .EQ. -999. .AND. &
k .GT. k_start .AND. qg (i,k,j) .GT. 1.E-20 ) THEN
w_n15 = w_2 (i,k,j)
qg_n15 = 1000. * qg (i,k,j) ! g/kg
grpl_flx = qg_n15 * w_n15
ENDIF
ENDDO
IF ( grpl_flx .GT. grpl_flx_max(i,j) ) THEN
grpl_flx_max(i,j) = grpl_flx
ENDIF
ENDDO
ENDDO
! Update buoyancy parameters.
! ---------------------------
IF ( do_buoy_calc ) THEN
nz = k_end - k_start + 1
DO j = jts, jte
DO i = its, ite
ostat = Buoyancy ( nz &
, t_phy(i,kms:kme ,j) &
, rh(i,kms:kme ,j) &
, p(i,kms:kme ,j) &
, z(i,kms:kme ,j) &
, 1 &
, cape(i,j) &
, cin(i,j) &
, zlfc_msl &
, plfc(i,j) &
, lidx(i,j) &
, 3 ) !Surface
! Add 273.15 to lifted index per some standard I don't understand
! ---------------------------------------------------------------
IF ( lidx ( i, j ) .ne. 999. ) lidx ( i, j ) = lidx ( i, j ) + 273.15
! Subtract terrain height to convert ZLFC from MSL to AGL
! -------------------------------------------------------
IF ( zlfc_msl .ge. 0. ) THEN
zlfc ( i, j ) = zlfc_msl - ht ( i, j )
ELSE
zlfc( i, j ) = -1.
ENDIF
ENDDO
ENDDO
ENDIF
! Maximum tornado wind speed in ms-1.
! First, save off old tornado mask and duration arrays
! ----------------------------------------------------
tornado_dur_prev(:,:)=tornado_dur(:,:)
tornado_mask_prev(:,:)=tornado_mask(:,:)
! Initialize all tornado variables
! --------------------------------
tornado_mask(:,:)=0.
tornado_dur(:,:)=0.
DO j = j_start, j_end
DO i = i_start, i_end
tornado = 0.
! Current tornado algorithm
! -------------------------
IF ( zlfc(i,j) .ge. 0. ) THEN
uh_term = min(max((uh(i,j) - 25.) / 50., 0.), 1.)
shr_term = min(max((llws(i,j) - 2.) / 10., 0.), 1.)
lfc_term = min(max((3000. - zlfc(i,j)) / 1500., 0.), 1.)
midrh2_term = min(max((90. - &
min(midrh_min(i,j),midrh_min_old(i,j))) / 30., 0.), 1.)
tornado = 50. * uh_term * shr_term * lfc_term * midrh2_term
ENDIF
! Does tornado algorithm indicate all possible ingredients
! for a minimum tornado, if so turn on mask
! --------------------------------------------------------
IF (tornado .GT. 0.) THEN
tornado_mask(i,j) = 1.
ENDIF
! Update the duration of this tornado if there was previously
! a tornado mask at this or an adjacent point
! -----------------------------------------------------------
IF ( ( tornado_mask(i,j) .GT. 0.5 ) .OR. &
( MAXVAL(tornado_mask_prev(i-1:i+1,j-1:j+1)) .GT. 0.5 ) ) THEN
tornado_dur(i,j) = MAXVAL(tornado_dur_prev(i-1:i+1,j-1:j+1)) + dt
ENDIF
! Ramp up value of tornado beta value linearly in time until &
! it has been sustained 5 minutes
! ------------------------------------------------------------
time_factor = MIN(tornado_dur(i,j)/900.,1.)
tornado = tornado * time_factor
! OPTIONAL
! Adjust previous max tornado wind upward to account for longer
! duration - if after 5 minutes, no adjustment made as
! time_factor/time_factor_prev=1
! -------------------------------------------------------------
!time_factor_prev = MIN((tornado_dur(i,j) - dt)/900.,1.)
!IF ( ( time_factor_prev .GT. 0. ) .AND. &
! ( time_factor_prev .LT. 1. ) ) THEN
! afwa_tornado(i,j) = afwa_tornado(i,j) * time_factor/time_factor_prev
!ENDIF
! Now that we are comparing apples to apples, see if current tornado
! wind is higher than previous highest value for this point
! ------------------------------------------------------------------
IF ( tornado .GT. afwa_tornado(i,j) ) THEN
afwa_tornado(i,j) = tornado
ENDIF
ENDDO
ENDDO
END SUBROUTINE severe_wx_diagnostics
SUBROUTINE vert_int_liquid_diagnostics ( vil & 1
, radarvil &
, t_phy &
, qrain &
, qsnow &
, qgrpl &
, z_e &
, dz8w &
, rho &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: rho &
, qrain &
, qsnow &
, qgrpl &
, t_phy &
, z_e &
, dz8w
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: vil &
, radarvil
! Local
! -----
INTEGER :: i,j,k,ktime
! Calculate vertically integrated liquid water (though its mostly not
! "liquid" now is it?) [kg/m^2]
! -------------------------------------------------------------------
DO i = ips, MIN(ipe,ide-1)
DO j = jps, MIN(jpe,jde-1)
vil (i,j) = 0.0
DO k = kps, MIN(kpe,kde-1)
vil (i,j) = vil (i,j) + &
(qrain (i,k,j) + &
qsnow (i,k,j) + &
qgrpl (i,k,j)) &
*dz8w (i,k,j) * rho(i,k,j)
ENDDO
ENDDO
ENDDO
! Diagnose "radar-derived VIL" from equivalent radar reflectivity
! radarVIL = (integral of LW*dz )/1000.0 (in kg/m^2)
! LW = 0.00344 * z_e** (4/7) in g/m^3
! ---------------------------------------------------------------
DO i = ips, MIN(ipe,ide-1)
DO j = jps, MIN(jpe,jde-1)
radarvil (i,j) = 0.0
DO k = kps, MIN(kpe,kde-1)
radarvil (i,j) = radarvil (i,j) + &
0.00344 * z_e(i,k,j)**0.57143 &
*dz8w (i,k,j)/1000.0
END DO
END DO
END DO
END SUBROUTINE vert_int_liquid_diagnostics
SUBROUTINE icing_diagnostics ( icing_opt & 1
, fzlev &
, icing_lg &
, icing_sm &
, qicing_lg_max &
, qicing_sm_max &
, qicing_lg &
, qicing_sm &
, icingtop &
, icingbot &
, t_phy &
, z &
, dz8w &
, rho &
, qrain &
, qcloud &
, ncloud &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
INTEGER, INTENT(IN) :: icing_opt
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: z &
, qrain &
, qcloud &
, ncloud &
, rho &
, dz8w &
, t_phy
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: fzlev &
, icing_lg &
, icing_sm &
, qicing_lg_max &
, qicing_sm_max &
, icingtop &
, icingbot
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT( OUT) :: qicing_lg &
, qicing_sm
! Local
! -----
INTEGER :: i,j,k,ktime,ktop,kbot
REAL :: qcfrac_lg, qcfrac_sm, qc, qr, small, all
! Initializations
! ---------------
fzlev (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
icingtop (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
icingbot (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
icing_lg (ips:ipe,jps:jpe) = 0.
icing_sm (ips:ipe,jps:jpe) = 0.
qicing_lg_max (ips:ipe,jps:jpe) = 0.
qicing_sm_max (ips:ipe,jps:jpe) = 0.
qicing_sm(ips:ipe,kps:kpe,jps:jpe)=0.
qicing_lg(ips:ipe,kps:kpe,jps:jpe)=0.
! Loop through i and j
! --------------------
DO i = ips, MIN(ipe,ide-1)
DO j = jps, MIN(jpe,jde-1)
! Go up the column and look for sub freezing temperatures
! -------------------------------------------------------
ktop=-1
kbot=-1
DO k = kps, MIN(kpe,kde-1)
IF (t_phy(i,k,j) .lt. 273.15) THEN
! Any cloud water we find will be supercooled.
! Based on microphysics scheme, determine the fraction of
! large (>50 um) supercooled cloud water drops.
! Source: Becky Selin, 16WS
! -------------------------------------------------------
qc = qcloud (i,k,j)
qr = qrain (i,k,j)
nc = ncloud(i,k,j)
den = rho(i,k,j)
qcfrac_lg = 0.
qcfrac_sm = 0.
! Eta (Ferrier)
! -------------
IF (icing_opt .eq. 2) THEN
IF (qc .lt. 2.5E-4) THEN
qcfrac_lg = 395000. * qc**2. + 102.9 * qc
ELSEIF (qc .lt. 1.4E-3) THEN
qcfrac_lg = 276.1 * qc - 0.01861
ELSE
qcfrac_lg = 0.3 * log(641.789 * qc) + 0.4
ENDIF
! Thompson
! --------
! RAS Per James McCormick's stats, more large supercooled
! drops are needed from the Thompson members. Changing
! calculation to be like WSM5/6 members.
!ELSEIF (icing_opt .eq. 3) THEN
! IF (qc .lt. 1.0E-3) THEN
! qcfrac_lg = 2205.0 * qc**2. + 3.232 * qc
! ELSEIF (qc .lt. 3.0E-3) THEN
! qcfrac_lg = 24.1 * qc - 0.01866
! ELSE
! qcfrac_lg = 0.127063 * log(550.0 * qc) - 0.01
! ENDIF
! Thompson or WSM5/6
! ------------------
!ELSEIF (icing_opt .eq. 4) THEN
ELSEIF ((icing_opt .eq. 3) .OR. (icing_opt .eq. 4)) THEN
IF (qc .lt. 5.E-4) THEN
qcfrac_lg = 50420.0 * qc**2. + 29.39 * qc
ELSEIF (qc .lt. 1.4E-3) THEN
qcfrac_lg = 97.65 * qc - 0.02152
ELSE
qcfrac_lg = 0.2 * log(646.908 * qc) + 0.135
ENDIF
! Morrison 2-moment, constant CCN
! -------------------------------
ELSEIF (icing_opt .eq. 5) THEN
IF (qc .lt. 1.4E-3) THEN
qcfrac_lg = 28000. * qc**2. + 0.1 * qc
ELSEIF (qc .lt. 2.6E-3) THEN
qcfrac_lg = 112.351 * qc - 0.102272
ELSE
qcfrac_lg = 0.3 * log(654.92 * qc) * 0.301607
ENDIF
! WDM6 or Morrison 2-moment w/ prognostic CCN
! -------------------------------------------
ELSEIF ((icing_opt .eq. 6) .OR. (icing_opt .eq. 7)) THEN
IF ((qc .gt. 1.0E-12) .and. (nc .gt. 1.0E-12)) THEN
small = -nc * exp(-nc*3141.59265*(5.E-5)**3./(6000.*den*qc))+nc
all = -nc * exp(-nc*3141.59265*(2.)**3./(6000.*den*qc))+nc
qcfrac_lg = 1. - (small / all)
ELSE
qcfac_lg = 0.
ENDIF
ENDIF
qcfrac_lg = max(qcfrac_lg, 0.)
! Small (<50 um) supercooled cloud water drop fraction (1 - large).
! -----------------------------------------------------------------
IF (icing_opt .ne. 0 ) THEN
qcfrac_sm = 1 - qcfrac_lg
ENDIF
! Supercooled drop mixing ratio
! -----------------------------
qicing_lg (i,k,j) = max(qr + qcfrac_lg * qc, 0.)
qicing_sm (i,k,j) = max(qcfrac_sm * qc, 0.)
! Column integrated icing
! -----------------------
icing_lg (i,j) = icing_lg (i,j) + qicing_lg (i,k,j) &
* dz8w (i,k,j) * rho(i,k,j)
icing_sm (i,j) = icing_sm (i,j) + qicing_sm (i,k,j) &
* dz8w (i,k,j) * rho(i,k,j)
! Column maximum supercooled drop mixing ratio
! --------------------------------------------
IF ( qicing_lg(i,k,j) .gt. qicing_lg_max(i,j) ) THEN
qicing_lg_max (i,j) = qicing_lg(i,k,j)
ENDIF
IF ( qicing_sm(i,k,j) .gt. qicing_sm_max(i,j) ) THEN
qicing_sm_max (i,j) = qicing_sm(i,k,j)
ENDIF
! Freezing level calculation
! --------------------------
IF (fzlev (i,j) .eq. -999.) THEN ! At freezing level
IF (k .ne. kps) THEN ! If not at surface, interpolate.
fzlev (i,j) = z (i,k-1,j) + &
((273.15 - t_phy (i,k-1,j)) &
/(t_phy (i,k,j) - t_phy (i,k-1,j))) &
*(z (i,k,j) - z (i,k-1,j))
ELSE ! If at surface, use first level.
fzlev(i,j) = z (i,k,j)
ENDIF
ENDIF
! Icing layer top and bottom indices (where icing > some arbitrary
! small value). Set bottom index of icing layer to current k index
! if not yet set. Set top index of icing layer to current k index.
! ----------------------------------------------------------------
IF ((qicing_lg (i,k,j) + qicing_sm (i,k,j)) .ge. 1.E-5) THEN
IF (kbot .eq. -1) kbot = k
ktop=k
ENDIF
ENDIF
END DO
! Interpolate bottom of icing layer from kbot (bottom index of icing
! layer). Icing bottom should not go below freezing level.
! ------------------------------------------------------------------
IF (kbot .ne. -1) THEN
IF (kbot .ne. kps) THEN ! If not at surface, interpolate
icingbot (i,j) = z (i,kbot-1,j) + ((1.E-5 - &
(qicing_lg (i,kbot-1,j) + qicing_sm (i,kbot-1,j))) &
/ ((qicing_lg (i,kbot,j) + qicing_sm (i,kbot,j)) &
- (qicing_lg (i,kbot-1,j) + qicing_sm (i,kbot-1,j)))) &
* (z (i,kbot,j) - z (i,kbot-1,j))
icingbot (i,j) = MAX(icingbot (i,j), fzlev (i,j))
ELSE ! If at surface use first level.
icingbot (i,j) = z(i,kbot,j)
ENDIF
ENDIF
! Interpolate top of icing layer from ktop (top index of icing layer).
! Icing top should not go below icing bottom (obviously).
! --------------------------------------------------------------------
IF (ktop .ne. -1 .and. ktop .ne. kpe) THEN ! If not undefined or model top
icingtop (i,j) = z (i,ktop,j) + ((1.E-5 - &
(qicing_lg (i,ktop,j) + qicing_sm (i,ktop,j))) &
/ ((qicing_lg (i,ktop+1,j) + qicing_sm (i,ktop+1,j)) &
- (qicing_lg (i,ktop,j) + qicing_sm (i,ktop,j)))) &
* (z (i,ktop+1,j) - z (i,ktop,j))
icingtop (i,j) = MAX(icingtop (i,j), icingbot (i,j))
ENDIF
END DO
END DO
END SUBROUTINE icing_diagnostics
SUBROUTINE radar_diagnostics ( refd & 1
, refd_com &
! , refd_max &
, echotop &
, z &
, z_e &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: z &
, z_e
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: refd &
, refd_com &
! , refd_max &
, echotop
! Local
! -----
INTEGER :: i,j,k,ktime
DO j = jps, MIN(jpe,jde-1)
DO i = ips, MIN(ipe,ide-1)
ktop = -1 ! Undefined
echotop (i,j) = 0.
refd_com (i,j) = 0.
refd (i,j) = 0.
DO k = kps, MIN(kpe,kde-1)
IF (z_e(i,k,j) .gt. 1.e-20) THEN
! Reflectivity (first level)
! --------------------------
IF (k == kps) refd(i,j) = MAX(10.0 * log10(z_e(i,k,j)),0.)
! ! Max reflectivity over the output interval
! ! -----------------------------------------
! IF (refd(i,j) .gt. refd_max(i,j)) refd_max(i,j) = refd(i,j)
! Composite reflectivity calc (max reflectivity in the column)
! ------------------------------------------------------------
IF (10.0 * log10(z_e(i,k,j)) .gt. refd_com(i,j)) THEN
refd_com(i,j) = 10.0 * log10(z_e(i,k,j))
ENDIF
ENDIF
! Echo top - the highest level w/ dBZ > 18 (z_e > 63.0957)
! --------------------------------------------------------
IF ( z_e(i,k,j) .gt. 63.0957) THEN
ktop = k
ENDIF
END DO
IF ( ktop .ne. -1 ) THEN ! Interpolate to echo top height (GAC)
echotop (i,j) = z (i,ktop,j) + &
((63.0957 - z_e (i,ktop,j)) &
/(z_e (i,ktop+1,j) - z_e (i,ktop,j))) &
*(z (i,ktop+1,j) - z (i,ktop,j))
ENDIF
END DO
END DO
END SUBROUTINE radar_diagnostics
SUBROUTINE wrf_dbzcalc( rho & 1
, t_phy &
, qr &
, qs &
, qg &
, z_e &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: rho &
, t_phy &
, qr &
, qs &
, qg
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT( OUT) :: z_e
REAL :: factor_r, factor_s, factor_g, factorb_s, factorb_g, ronv, sonv, gonv
REAL :: temp_c, rhoair, qgr, qra, qsn
INTEGER :: i, j, k
INTEGER, PARAMETER :: iBrightBand = 1
REAL, PARAMETER :: T_0 = 273.15
REAL, PARAMETER :: PI = 3.1415926536
REAL, PARAMETER :: rgas=287.04, gamma_seven = 720.0, alpha2 = 0.224
! Densities of rain, snow, graupel, and cloud ice.
! ------------------------------------------------
REAL, PARAMETER :: rho_w = 1000.0, rho_r = 1000.0, rho_s = 100.0
REAL, PARAMETER :: rho_g = 400.0, rho_i = 890.0
REAL, PARAMETER :: ron=8.e6, son=2.e7, gon=5.e7, r1=1.e-15
REAL, PARAMETER :: ron_min = 8.e6, ron2=1.e10
REAL, PARAMETER :: ron_qr0 = 0.0001, ron_delqr0 = 0.25*ron_qr0
REAL, PARAMETER :: ron_const1r = (ron2-ron_min)*0.5
REAL, PARAMETER :: ron_const2r = (ron2+ron_min)*0.5
! Constant intercepts
! -------------------
ronv = 8.e6 ! m^-4
sonv = 2.e7 ! m^-4
gonv = 4.e6 ! m^-4
factor_r = gamma_seven * 1.e18 * (1./(pi*rho_r))**1.75
factor_s = gamma_seven * 1.e18 * (1./(pi*rho_s))**1.75 &
* (rho_s/rho_w)**2 * alpha2
factor_g = gamma_seven * 1.e18 * (1./(pi*rho_g))**1.75 &
* (rho_g/rho_w)**2 * alpha2
! For each grid point
! -------------------
DO j = jps, jpe
DO k = kps, kpe
DO i = ips, ipe
factorb_s = factor_s
factorb_g = factor_g
! In this case snow or graupel particle scatters like liquid
! water because it is assumed to have a liquid skin
! ----------------------------------------------------------
IF( iBrightBand == 1 ) THEN
IF (t_phy(i,k,j) > T_0) THEN
factorb_s = factor_s /alpha2
factorb_g = factor_g /alpha2
ENDIF
ENDIF
! Calculate variable intercept parameters
! ---------------------------------------
temp_c = amin1(-0.001, t_phy(i,k,j)- T_0)
sonv = amin1(2.0e8, 2.0e6*exp(-0.12*temp_c))
gonv = gon
qgr = QG(i,k,j)
qra = QR(i,k,j)
qsn = QS(i,k,j)
IF (qgr.gt.r1) THEN
gonv = 2.38*(pi*rho_g/(rho(i,k,j)*qgr))**0.92
gonv = max(1.e4, min(gonv,gon))
ENDIF
ronv = ron2
IF (qra.gt. r1) THEN
ronv = ron_const1r*tanh((ron_qr0-qra)/ron_delqr0) + ron_const2r
ENDIF
IF (qra < 0.0 ) qra = 0.0
IF (qsn < 0.0 ) qsn = 0.0
IF (qgr < 0.0 ) qgr = 0.0
z_e(i,k,j) = factor_r * (rho(i,k,j) * qra)**1.75 / ronv**.75 + &
factorb_s * (rho(i,k,j) * qsn)**1.75 / sonv**.75 + &
factorb_g * (rho(i,k,j) * qgr)**1.75 / gonv**.75
IF ( z_e(i,k,j) < 0.0 ) z_e(i,k,j) = 0.0
END DO
END DO
END DO
END SUBROUTINE wrf_dbzcalc
SUBROUTINE precip_type_diagnostics ( t_phy & 1
, t2 &
, rh &
, z &
, ht &
, precip &
, swdown &
, rain &
, snow &
, ice &
, frz_rain &
, snowfall &
, ccn_tmp &
, total_melt &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: t_phy &
, rh &
, z
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: t2 &
, ht &
, precip &
, swdown
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT) :: snowfall &
, rain &
, frz_rain &
, snow &
, ice
REAL, INTENT(IN) :: ccn_tmp
REAL, INTENT(IN) :: total_melt
! Local
! -----
REAL, DIMENSION( ims:ime, jms:jme ) :: &
melt &
, mod_2m_tmp &
, cloud_top_tmp &
, maxtmp
INTEGER, DIMENSION( ims:ime, jms:jme ) :: &
cloud_top_k_index &
, precip_type
LOGICAL, DIMENSION (ims:ime, jms:jme ) :: &
saturation
REAL, PARAMETER :: snow_ratio=5.0
! Loop through all points
! Search vertically twice--first to find the cloud top temperature and the
! maximum temperature. Second, determine if any melting or re-freezing will
! occur to make ice pellets or freezing rain
! -------------------------------------------------------------------------
DO i=ips,ipe
DO j=jps,jpe
saturation(i,j)=.false.
melt(i,j)=0.0
precip_type(i,j)=0
! Modify surface temperature for solar insolation (W/m2)
! Set max temperature in the atmopshere
! ------------------------------------------------------
mod_2m_tmp(i,j)=t2(i,j)+(swdown(i,j)/100.0)
maxtmp(i,j)=mod_2m_tmp(i,j)
! Only look at points that have precip and are not warm at the surface
! --------------------------------------------------------------------
IF (precip(i,j) .gt. 0.0) THEN
!IF (mod_2m_tmp(i,j) .gt. 277.15) THEN
IF (mod_2m_tmp(i,j) .gt. 275.15) THEN
precip_type(i,j)=1 ! Rain
ELSE
! Check sounding from top for saturation (RH-water gt 80%)--this is
! the cloud top. Erase saturation if RH lt 70% (spurious moist layer
! aloft)
! ------------------------------------------------------------------
cloud_top_k_index(i,j)=kpe
DO k=kpe,kps,-1
IF ((z(i,k,j)-ht(i,j)) .gt. 0.0) THEN
IF (t_phy(i,k,j) .gt. maxtmp(i,j)) THEN
maxtmp(i,j)=t_phy(i,k,j)
ENDIF
IF (rh(i,k,j) .gt. 80 .and. saturation(i,j) .eqv. .false.) THEN
cloud_top_tmp(i,j)=t_phy(i,k,j)
cloud_top_k_index(i,j)=k
saturation(i,j)=.true.
precip_type(i,j)=2 ! Snow
ENDIF
IF (rh(i,k,j) .le. 70 .and. saturation(i,j) .eqv. .true.) THEN
saturation(i,j)=.false.
ENDIF
ENDIF
ENDDO
! Perform simple check to assign types with no melting layer
! shenanigans going on
! ----------------------------------------------------------------
IF (cloud_top_tmp(i,j) .le. ccn_tmp .and. &
maxtmp(i,j) .le. 273.15) THEN
precip_type(i,j)=2 ! Snow
ENDIF
! ELSE, have to go through the profile again to see if snow melts,
! and if anything re-freezes
! ----------------------------------------------------------------
DO k=cloud_top_k_index(i,j),kps,-1
IF ((z(i,k,j)-ht(i,j)) .gt. 0.0) THEN
! Condition 0--assign falling rain when we get to the
! supercooled temperature if too warm
! ---------------------------------------------------
IF (cloud_top_tmp(i,j) .eq. t_phy(i,k,j) .and. &
cloud_top_tmp(i,j) .gt. ccn_tmp) THEN
precip_type(i,j)=1 ! Rain
ENDIF
! Condition 1--falling frozen precip that will start to melt
! Add up melting energy over warm layers--if enough, turn to
! liquid
! ----------------------------------------------------------
IF ((precip_type(i,j) .eq. 2 .or. precip_type(i,j) .eq. 3) .and. &
t_phy(i,k,j) .gt. 273.15) THEN
melt(i,j)=melt(i,j)+9.8*(((t_phy(i,k,j)-273.15)/273.15)* &
(z(i,k,j)-z(i,k-1,j)))
IF (melt(i,j) .gt. total_melt) THEN
precip_type(i,j)=1 ! Rain
melt(i,j)=0.0 ! Reset melting energy in case it re-freezes
ENDIF
ENDIF
! Condition 2--falling partially melted precip encounters
! sub-freezing air. Snow will be converted to ice pellets if
! at least 1/4 of it melted. Instantaneous freeze-up, simplistic
! --------------------------------------------------------------
IF (t_phy(i,k,j) .le. 273.15 .and. &
melt(i,j) .gt. total_melt/4.0 .and. &
(precip_type(i,j) .eq. 2 .or. precip_type(i,j) .eq. 3)) THEN
precip_type(i,j)=3 ! Ice
melt(i,j)=0.0
ENDIF
! Condition 3--falling liquid that will re-freeze--must reach
! nucleation temperature
! -----------------------------------------------------------
IF (precip_type(i,j) .eq. 1) THEN
IF (t_phy(i,k,j) .le. ccn_tmp) THEN
precip_type(i,j)=3 ! Ice
ENDIF
ENDIF
ENDIF ! End if (z-ht)>0
ENDDO ! End do k=kpe,kps,-1
ENDIF ! End if mod_2m_tmp>273.15
! Accumulate precip according to precip_type
! ------------------------------------------
IF (precip_type(i,j) .eq. 3) THEN
ice(i,j)=ice(i,j)+precip(i,j)
ENDIF
IF (precip_type(i,j) .eq. 2) THEN
snow(i,j)=snow(i,j)+precip(i,j)
snowfall(i,j)=snowfall(i,j)+snow_ratio*precip(i,j) &
*(5.-mod_2m_tmp(i,j)+273.15)**0.4
ENDIF
IF (precip_type(i,j) .eq. 1) THEN
IF (mod_2m_tmp(i,j) .gt. 273.15) THEN
rain(i,j)=rain(i,j)+precip(i,j)
ELSE
frz_rain(i,j)=frz_rain(i,j)+precip(i,j)
ENDIF
ENDIF
ENDIF ! End if precip>0
ENDDO ! End do j=jps,jpe
ENDDO ! End do i=ips,ipe
END SUBROUTINE precip_type_diagnostics
SUBROUTINE vis_diagnostics ( qcloud & 1
, qrain &
, qice &
, qsnow &
, qgrpl &
, rho &
, wind10m &
, wind125m &
, pwater &
, q2m &
, rh2m &
, rh20m &
, tv2m &
, tv20m &
, dustc &
, vis &
, vis_dust &
, vis_alpha &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
INTEGER, PARAMETER :: ndust=5
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: qcloud &
, qrain &
, qice &
, qsnow &
, qgrpl &
, rho &
, wind10m &
, wind125m &
, pwater &
, rh2m &
, q2m &
, rh20m &
, tv2m &
, tv20m
REAL, DIMENSION( ims:ime, jms:jme, ndust ), &
INTENT(IN ) :: dustc
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: vis &
, vis_dust &
, vis_alpha
! Local
! -----
INTEGER :: i,j,k,d
REAL, PARAMETER :: visfactor=3.912
REAL, DIMENSION (ndust) :: dustfact
REAL :: bc, br, bi, bs, dust_extcoeff, hydro_extcoeff, extcoeff, vis_haze
REAL :: tvd, rh, prob_ext_coeff_gt_p29, haze_ext_coeff
REAL :: vis_hydlith, alpha_haze
! Dust factor based on 5 bin AFWA dust scheme. This is a simplification
! of the scheme in WRFPOST. More weight is applied to smaller particles.
! -----------------------------------------------------------------------
dustfact=(/1.470E-6,7.877E-7,4.623E-7,2.429E-7,1.387E-7/)
DO i=ims,ime
DO j=jms,jme
! Hydrometeor extinction coefficient
! For now, lump graupel in with rain
! -------------------------------------------------------------------
! Update: GAC 20131213 Our test results with surface cloud and ice
! are very unfavorable. Model doesn't do a great job handling clouds
! at the model surface. Therefore, we will not trust surface
! cloud water/ice. (Commented out below).
! -------------------------------------------------------------------
!br=2.240*qrain(i,j)**0.75
!bs=10.36*qsnow(i,j)**0.78
!bc=144.7*qcloud(i,j)**0.88
!bi=327.8*qice(i,j)
!hydro_extcoeff=bc+br+bi+bs
!br=2.240*(qrain(i,j)+qgrpl(i,j))**0.75
!bs=10.36*(qsnow(i,j)*rho(i,j))**0.78
! Update: moisture variables should be in mass concentration (g m^-3)
br=1.1*(1000.*rho(i,j)*(qrain(i,j)+qgrpl(i,j)))**0.75
bs=10.36*(1000.*rho(i,j)*qsnow(i,j))**0.78
hydro_extcoeff=(br+bs)/1000. ! m^-1
! Dust extinction coefficient
! ---------------------------
dust_extcoeff=0.
DO d=1,ndust
dust_extcoeff=dust_extcoeff+dustfact(d)*dustc(i,j,d)
ENDDO
! UPDATE: GAC 20131213 Old algorithm commented out below
!! Visibility due to haze obscuration
!! -------------------------------------------------------
!vis_haze=1500.*(105.-rh2m(i,j)+wind10m(i,j))
!
!! Calculate total visibility
!! Take minimum visibility from hydro/lithometeors and haze
!! Define maximum visibility as 20 km (UPDATE: 999.999 km)
!! --------------------------------------------------------
!extcoeff=hydro_extcoeff+dust_extcoeff
!IF (extcoeff .gt. 0.) THEN
! vis(i,j)=MIN(visfactor/extcoeff,vis_haze)
!ELSE
! vis(i,j)=999999.
!ENDIF
! Update: GAC 20131213 New haze/fog visibility algorithm
! Start with relative humidity predictor. Increase
! predicted visibility as mixing ratio decreases (as
! there is less water to condense).
! -------------------------------------------------------
vis_haze=999999.
IF (q2m(i,j) .gt. 0.) THEN
!vis_haze=1500.*(105.-rh2m(i,j))*(15./(1000.*q2m(i,j)))
vis_haze=1500.*(105.-rh2m(i,j))*(5./min(1000.*q2m(i,j),5.))
ENDIF
! Calculate a Weibull function "alpha" term. This can be
! used later with visibility (which acts as the "beta" term
! in the Weibull function) to create a probability distribution
! for visibility. Alpha can be thought of as the "level of
! certainty" that beta (model visibility) is correct. Fog is
! notoriously difficult to model. In the below algorithm,
! the alpha value (certainty) decreases as PWAT, mixing ratio,
! or winds decrease (possibly foggy conditions), but increases
! if RH decreases (more certainly not foggy). If PWAT is lower
! then there is a higher chance of radiation fog because there
! is less insulating cloud above.
! -------------------------------------------------------------
alpha_haze=3.6
IF (q2m(i,j) .gt. 0.) THEN
alpha_haze=0.1 + pwater(i,j)/25. + wind125m(i,j)/3. + &
(100.-rh2m(i,j))/10. + 1./(1000.*q2m(i,j))
alpha_haze=min(alpha_haze,3.6)
ENDIF
! Calculate visibility from hydro/lithometeors
! Maximum visibility -> 999999 meters
! --------------------------------------------
extcoeff=hydro_extcoeff+dust_extcoeff
IF (extcoeff .gt. 0.) THEN
vis_hydlith=min(visfactor/extcoeff, 999999.)
ELSE
vis_hydlith=999999.
ENDIF
! Calculate total visibility
! Take minimum visibility from hydro/lithometeors and haze
! Set alpha to be alpha_haze if haze dominates, or 3.6
! (a Guassian distribution) when hydro/lithometeors dominate
! ----------------------------------------------------------
IF (vis_hydlith < vis_haze) THEN
vis(i,j)=vis_hydlith
vis_alpha(i,j)=3.6
ELSE
vis(i,j)=vis_haze
vis_alpha(i,j)=alpha_haze
ENDIF
! Calculate dust visibility
! Again, define maximum visibility as 20 km
! -----------------------------------------
IF (dust_extcoeff .gt. 0.) THEN
vis_dust(i,j)=MIN(visfactor/dust_extcoeff,999999.)
ELSE
vis_dust(i,j)=999999.
ENDIF
ENDDO
ENDDO
END SUBROUTINE vis_diagnostics
SUBROUTINE cloud_diagnostics (qcloud & 1
, qice &
, qsnow &
, rh &
, dz8w &
, rho &
, z &
, ht &
, cloud &
, cloud_ceil &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: qcloud &
, qice &
, qsnow &
, rh &
, dz8w &
, rho &
, z
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: ht
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: cloud &
, cloud_ceil
! Local
! -----
INTEGER :: i, j, k
REAL :: tot_cld_cond, maxrh, cld_frm_cnd, cld_frm_rh, z_maxrh
REAL :: snow_extcoeff, vis_snow, cloud_lo, zagl_up, zagl_lo
REAL, PARAMETER :: min_ceil = 125. ! Minimum ceiling of 125 m
! Calculate cloud cover based on total cloud condensate, or if none
! present, from maximum relative humidity in the column.
! -----------------------------------------------------------------
DO i=ims,ime
DO j=jms,jme
! Initialize some key variables
! -----------------------------
tot_cld_cond = 0.
cloud(i,j) = 0.
maxrh = -9999.
cloud_ceil(i,j) = -9999.
cloud_lo = 0.
! Now go up the column to find our cloud base
! -------------------------------------------
DO k=kms,kme
! Total cloud condensate
! Don't trust modeled cloud below 125 m. We
! let the alpha curve determine probabilities
! below 125 m...
! ---------------------------------------------
IF ( z(i,k,j) - ht (i,j) .gt. min_ceil ) THEN
! Maximum column relative humidity
! --------------------------------
IF (rh (i,k,j) .gt. maxrh) THEN
maxrh = rh (i,k,j)
z_maxrh = z(i,k,j)
ENDIF
! ! Cloud cover parameterization. Take maximum value
! ! from condensate and relative humidity terms.
! ! ------------------------------------------------
! tot_cld_cond = tot_cld_cond + (qcloud (i,k,j) + qice (i,k,j) &
! + qsnow (i,k,j)) * dz8w (i,k,j) * rho(i,k,j)
! cld_frm_cnd = 50. * tot_cld_cond
! cld_frm_rh = MAX(((maxrh - 70.) / 30.),0.)
! cloud (i,j) = MAX(cld_frm_cnd,cld_frm_rh)
! Calculate cloud cover beta value by summing
! relative humidity above 70% as we go up the
! column. Assume a higher probability of a
! cloud if we have an accumulation of high
! relative humidity over a typical cloud
! depth of 500m. (Note dz8w is distance
! between full eta levels. Note also that rh
! is derived from the sum of qvapor, qcloud,
! and qcloud, with a maximum of 100%).
! -------------------------------------------
cld_frm_rh = MAX(((rh (i,k,j) - 90.) / 10.),0.)
cloud (i,j) = cloud (i,j) + ( cld_frm_rh * dz8w (i,k,j) ) / 250.
! Calculate cloud ceiling, the level at which
! parameterization of cloud cover exceeds 80%
! Once we exceed the 80% threshold, we will
! interpolate downward to find the ceiling.
! If this is the lowest level, then we will
! simply set ceiling to that level. After
! we interpolate, if ceiling is below the
! minimum we trust, we set it to the minimum.
! -------------------------------------------
IF ( cloud_ceil (i,j) .eq. -9999. .and. cloud (i,j) .gt. 0.8 ) THEN
zagl_up = z (i,k,j) - ht (i,j)
IF ( k .EQ. kps ) THEN
cloud_ceil (i,j) = zagl_up
ELSE
zagl_lo = z (i,k-1,j) - ht (i,j)
cloud_ceil (i,j) = zagl_lo + &
((0.8 - cloud_lo) / &
(cloud (i,j) - cloud_lo)) * &
(zagl_up - zagl_lo)
cloud_ceil (i,j) = MAX(cloud_ceil (i,j),ceil_min)
ENDIF
ENDIF
! Save cloud amount here to use for interpolation
! -----------------------------------------------
cloud_lo=cloud(i,j)
ENDIF
ENDDO
! If we did not encounter any definitive cloud earlier, we set
! cloud ceiling to the level of maximum relative humidity. (If
! there is any cloud, this is our best guess as to where it
! will reside in the vertical). Height is AGL.
! ------------------------------------------------------------
IF (cloud_ceil (i,j) .eq. -9999.) THEN
cloud_ceil (i,j) = z_maxrh - ht (i,j)
ENDIF
! Compare cloud ceiling to surface visibility reduction due to snow
! Note difference from horizontal visibility algorithm for snow
! -----------------------------------------------------------------
IF (qsnow (i,1,j) .GT. 0. .AND. rho (i,1,j) .GT. 0.) THEN
snow_extcoeff = 25. * (1000. * rho(i,1,j) * qsnow (i,1,j))**0.78
snow_extcoeff = snow_extcoeff / 1000.
vis_snow = 3.912 / snow_extcoeff
IF (vis_snow .LT. cloud_ceil (i,j)) cloud_ceil (i,j) = vis_snow
ENDIF
ENDDO
ENDDO
END SUBROUTINE cloud_diagnostics
SUBROUTINE thermal_diagnostics ( t2 & 1
, psfc &
, rh2m &
, wind10m &
, heatidx &
, wchill &
, fits &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: t2 &
, psfc &
, rh2m &
, wind10m
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: heatidx &
, wchill &
, fits
! Local
! -----
INTEGER :: i, j
DO i=ims,ime
DO j=jms,jme
! Heat Index
! ----------
heatidx ( i, j ) = calc_hi ( t2 ( i, j ) &
, rh2m ( i, j ) )
! Wind Chill
! ----------
wchill ( i, j ) = calc_wc ( t2 ( i, j ) &
, wind10m ( i, j ) )
! Fighter Index of Thermal Stress
! -------------------------------
fits ( i, j ) = calc_fits ( psfc ( i, j ) &
, t2 ( i, j ) &
, rh2m ( i, j ) )
ENDDO
ENDDO
END SUBROUTINE thermal_diagnostics
SUBROUTINE turbulence_diagnostics ( u_phy & 1
, v_phy &
, t_phy &
, p &
, zagl &
, defor11 &
, defor12 &
, defor22 &
, turb &
, llturb &
, llturblgt &
, llturbmdt &
, llturbsvr &
, nlyrs &
, lyrbot &
, lyrtop &
, ids, ide, jds, jde, kds, kde &
, ims, ime, jms, jme, kms, kme &
, ips, ipe, jps, jpe, kps, kpe )
INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
INTEGER, INTENT(IN) :: nlyrs
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: u_phy &
, v_phy &
, t_phy &
, p &
, zagl &
, defor11 &
, defor12 &
, defor22
REAL, DIMENSION( nlyrs ), &
INTENT(IN ) :: lyrtop &
, lyrbot
REAL, DIMENSION( ims:ime, nlyrs, jms:jme ), &
INTENT( OUT) :: turb
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT( OUT) :: llturb &
, llturblgt &
, llturbmdt &
, llturbsvr
! Local
! -----
INTEGER :: i, j, k, n, bot, top, nlayer
REAL :: ugrdtop &
, ugrdbot &
, vgrdtop &
, vgrdbot &
, defor11top &
, defor11bot &
, defor12top &
, defor12bot &
, defor22top &
, defor22bot
REAL, DIMENSION( kms:kme ) :: this_zagl &
, this_tK &
, this_p &
, this_u &
, this_v
REAL :: wind, therm, mtn_wave, tpd_wave
!~ Initialize variables.
! ----------------------
turb = REAL ( 0 )
llturb = REAL ( 0 )
llturblgt = REAL ( 0 )
llturbsvr = REAL ( 0 )
!~ Loop through the grid.
! ----------------------
DO i=ims,ime
DO j=jms,jme
!~ Loop through the turbulence layers
! ----------------------------------
DO n = 1, nlyrs
!~ Interpolate relevent variables to turbulence layers
! ---------------------------------------------------
ugrdtop = lin_interp ( zagl ( i, kms:kme-1, j ) &
, u_phy ( i, kms:kme-1, j ) &
, lyrtop ( n ) )
ugrdbot = lin_interp ( zagl ( i, kms:kme-1, j ) &
, u_phy ( i, kms:kme-1, j ) &
, lyrbot ( n ) )
vgrdtop = lin_interp ( zagl ( i, kms:kme-1, j ) &
, v_phy ( i, kms:kme-1, j ) &
, lyrtop ( n ) )
vgrdbot = lin_interp ( zagl ( i, kms:kme-1, j ) &
, v_phy ( i, kms:kme-1, j ) &
, lyrbot ( n ) )
defor11top = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor11 ( i, kms:kme-1, j ) &
, lyrtop ( n ) )
defor11bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor11 ( i, kms:kme-1, j ) &
, lyrbot ( n ) )
defor12top = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor12 ( i, kms:kme-1, j ) &
, lyrtop ( n ) )
defor12bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor12 ( i, kms:kme-1, j ) &
, lyrbot ( n ) )
defor22top = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor22 ( i, kms:kme-1, j ) &
, lyrtop ( n ) )
defor22bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
, defor22 ( i, kms:kme-1, j ) &
, lyrbot ( n ) )
!~ Compute Knapp-Ellrod clear air turbulence
! -----------------------------------------
turb ( i, n, j ) = CATTurbulence ( ugrdbot &
, ugrdtop &
, vgrdbot &
, vgrdtop &
, defor11bot &
, defor11top &
, defor12bot &
, defor12top &
, defor22bot &
, defor22top &
, lyrbot (n) &
, lyrtop (n) )
ENDDO
!~ Get top and bottom index of 0-1500 m AGL layer
! ----------------------------------------------
bot = kms
top = kms
DO k=kms+1,kme
IF ( zagl ( i, k, j ) .gt. 1500. ) THEN
top = k
EXIT
ENDIF
ENDDO
nlayer = top - bot + 1 !~ Number of layers from bottom to top
!~ Copy current column at this i,j point into working arrays
! ---------------------------------------------------------
this_zagl = zagl ( i, kms:kme, j )
this_tK = t_phy ( i, kms:kme, j )
this_p = p ( i, kms:kme, j )
this_u = u_phy ( i, kms:kme, j )
this_v = v_phy ( i, kms:kme, j )
!~ Interpolate req'd vars to the 1500 m level
!~ Overwrite the "top" index with these values
! -------------------------------------------
this_zagl ( top ) = 1500.
this_tK ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
, t_phy ( i, kms:kme-1, j ) &
, this_zagl (top) )
this_p ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
, p ( i, kms:kme-1, j ) &
, this_zagl (top) )
this_u ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
, u_phy ( i, kms:kme-1, j ) &
, this_zagl (top) )
this_v ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
, v_phy ( i, kms:kme-1, j ) &
, this_zagl (top) )
!~ Compute the low level turbulence index (from 0 - 1500 m AGL)
!~ There are four components to this index: a wind speed term,
!~ thermodynamic term, mountain wave term, and trapped wave term.
!~ These terms will utilize the working arrays we have just
!~ defined, using only the portion of the array valid in the
!~ 0-1500 m layer, e.g. this_u (bot:top).
!~ The algorithm itself was developed by:
!~
!~ Mr. James McCormick
!~ Aviation Hazards Team
!~ Air Force Weather Agency 16WS/WXN
!~ DSN: 271-1689 Comm: (402) 294-1689
!~ James.McCormick.ctr@offutt.af.mil
! -------------------------------------------------------------- !
!~ Step 1: Compute the wind speed term ~!
! -------------------------------------------------------------- !
wind = LLT_WindSpeed ( nlayer, this_u (bot:top) &
, this_v (bot:top) )
! -------------------------------------------------------------- !
!~ Step 2: Compute the thermodynamic term ~!
! -------------------------------------------------------------- !
therm = LLT_Thermodynamic ( nlayer, this_tK(bot:top) &
, this_zagl(bot:top) )
! -------------------------------------------------------------- !
!~ Step 3: Compute the mountain wave term ~!
! -------------------------------------------------------------- !
mtn_wave = LLT_MountainWave ( nlayer, terrain_dx, terrain_dy &
, this_u(bot:top), this_v(bot:top) &
, this_tK(bot:top), this_zagl(bot:top) )
! -------------------------------------------------------------- !
!~ Step 4: Compute the trapped wave term ~!
! -------------------------------------------------------------- !
tpd_wave = LLT_TrappedWave ( nlayer, this_u(bot:top) &
, this_v(bot:top), this_p(bot:top) )
! -------------------------------------------------------------- !
!~ Step 5: Combine the above and arrive at the turbulence index. ~!
! -------------------------------------------------------------- !
llturb ( i,j ) = 1.-((1.-wind)*(1.-therm)*(1.-mtn_wave)*(1.-tpd_wave))
!~ Compute probabilities of light, moderate, and severe LLT
! --------------------------------------------------------
llturblgt ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (30)) &
*2.5)*.01)**2)*0.75)*REAL(100))
IF ( llturb (i,j) < 0.3 ) llturblgt ( i,j ) = REAL ( 0 )
IF ( llturblgt (i,j) > REAL (90) ) llturblgt ( i,j ) = REAL ( 90 )
llturbmdt ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (35)) &
*2.22222)*.01)*0.75)**2)*88.88888)
IF ( llturb (i,j) < 0.35 ) llturbmdt ( i,j ) = REAL ( 0 )
IF ( llturbmdt (i,j) > REAL (70) ) llturbmdt ( i,j ) = REAL ( 70 )
llturbsvr ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (40)) &
*REAL(2))*.01)*0.5)**2)*REAL(100))
IF ( llturb (i,j) < 0.40 ) llturbsvr ( i,j ) = REAL ( 0 )
IF ( llturbsvr (i,j) > REAL (35) ) llturbsvr ( i,j ) = REAL ( 35 )
ENDDO
ENDDO
END SUBROUTINE turbulence_diagnostics
END MODULE module_diag_afwa
#endif