Skip to content

Commit

Permalink
Bug fix in buoyancy diagnostic function (wrf-model#1430)
Browse files Browse the repository at this point in the history 
Correct MU layer and PWAT, refine CAPE, CIN, LFC

TYPE: bug fix, enhancement

KEYWORDS: Most unstable layer, Precipitable water, CAPE, CIN, LFC

SOURCE: Zhixiao Zhang (University of Utah), Adam Varble (PNNL and University of Utah), and Katelyn Barber (PNNL)

DESCRIPTION OF CHANGES:
Problem:
The most unstable (MU) layer and precipitable water(PWAT) are inconsistent with the conventional definitions in 
meteorological communities. LFC, CAPE and CIN calculations do not work well while dealing with multiple inversion 
layers.

Solution:
1. Correct the MU layer definition from the max water vapor mixing ratio (w) level to the max equivalent potential 
temperature (Theta-e) level. Because it is the theta-e rather than mixing ratio conserved during the dry/moist 
adiabatic lifting processes.
2. Correct PWAT definition from column accumulated total mass of water vapor and cloud water to only water vapor. 
This is based on the AMS glossary https://glossary.ametsoc.org/wiki/Precipitable_water.
3. Calculate MUCAPE, MUCIN and LFC by refined definitions: MUCAPE and MUCIN are defined as the vertically 
integrated positive and negative buoyant energy between LCL and EL with the air parcel lifted from MU layer. LFC 
is defined as the highest separation level when buoyancy turns from negative to positive.

LIST OF MODIFIED FILES:
phys/module_diag_functions.F

TESTS CONDUCTED: 
1. The modified code has been complied successfully and applied on regional climate simulations for CACTI field 
campaign.

2. In elevated deep convection regimes when the highest water vapor mixing ratios are still at low levels but the 
most unstable parcel is at a higher altitude, the differences can be large. An example of such a sounding is 
attached. The original diagnostic would choose a parcel starting at the surface where the highest water vapor 
mixing ratio exists. That parcel has no CAPE and obviously a lot of CIN. The most unstable parcel with the 
highest theta-e (you can find this using the highest dashed green moist adiabat value) is just above the 700-mb 
level which is a parcel with CAPE and very little CIN.
![upperair DOE_M1_sonde 201811121200 skewT](https://user-images.githubusercontent.com/17436975/111506933-5ffeff80-8707-11eb-89ff-a24590d419c2.png)

3. Jenkins testing is all PASS.

OTHER NOTES: The modification is based on WRF 4.1.1. The relevant AFWA module, diagnostic driver and registry require to be consistently modified in the latest source code by maintainers.

RELEASE NOTE: A few updates to the diagnostic schemes were added. The most unstable (MU) layer and precipitable water(PWAT) were inconsistent with the conventional definitions used by the meteorological communities. LFC, CAPE and CIN calculations previously did not work well while dealing with multiple inversion layers.
  • Loading branch information
@zhixiaozhang @vlakshmanan-scala
zhixiaozhang authored and vlakshmanan-scala committed Apr 4, 2024
1 parent 259d88f commit 9662d57
Showing 1 changed file with 61 additions and 47 deletions.
108 changes: 61 additions & 47 deletions phys/module_diag_functions.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -388,7 +388,10 @@ FUNCTION Pwat ( nz, qv, qc, dz8w, rho )
! ---------------------------
Pwat=0
DO k = 1, nz
Pwat = Pwat + (qv(k) + qc(k)) * dz8w(k) * rho(k)
!Based on AMS PWAT defination (https://glossary.ametsoc.org/wiki/Precipitable_water)
!PWAT is corrected as the column accumulated water vapor rather than water vapor + cloud water.
!Modified by Zhixiao
Pwat = Pwat + qv(k) * dz8w(k) * rho(k)
ENDDO
END FUNCTION Pwat
Expand Down Expand Up @@ -463,10 +466,12 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
! -------------------------
REAL :: ws ( nz ) !~ Saturation mixing ratio
REAL :: w ( nz ) !~ Mixing ratio
REAL :: etheta( nz )!~ Equivalent potential temperature. Modified by Zhixiao.
REAL :: dTvK ( nz ) !~ Parcel / ambient Tv difference
REAL :: buoy ( nz ) !~ Buoyancy
REAL :: tlclK !~ LCL temperature ( K )
REAL :: plcl !~ LCL pressure ( Pa )
REAL :: pel !~ Equilibrium pressure ( Pa ). Modified by Zhixiao.
REAL :: nbuoy !~ Negative buoyancy
REAL :: pbuoy !~ Positive buoyancy

Expand All @@ -481,6 +486,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
REAL :: srcthetaeK !~ Source parcel theta-e ( K )
INTEGER :: srclev !~ Level of the source parcel
REAL :: spdiff !~ Pressure difference
REAL :: srce !~ Equivalent potential temperature ( K ). Modified by Zhixiao.

!~ Parcel variables
! ----------------
Expand All @@ -493,6 +499,8 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
! -----------------------
INTEGER :: i, j, k !~ Dummy iterator
INTEGER :: lfclev !~ Level of LFC
INTEGER :: ellev !~ Level of EL. Modified by Zhixiao.
INTEGER :: lcllev !~ Level of LCL. Modified by Zhixiao.
INTEGER :: prcl !~ Internal parcel type indicator
INTEGER :: mlev !~ Level for ML calculation
INTEGER :: lyrcnt !~ Number of layers in mean layer
Expand All @@ -514,12 +522,11 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
PARAMETER ( g = 9.80665 ) !~ m s^-2
REAL :: RUNDEF
PARAMETER ( RUNDEF = -9.999E30 )

!~ Initialize variables
! --------------------
ostat = 0
CAPE = REAL ( 0 )
CIN = REAL ( 0 )
CIN = RUNDEF !Change CIN filling values from 0 to default filling. CIN should not initially be filled by 0, because 0 means no inhibition energy. Modified by Zhixiao
ZLFC = RUNDEF
PLFC = RUNDEF

Expand Down Expand Up @@ -549,6 +556,11 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
DO k = sfc, nz
ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
w ( k ) = ( rh(k)/100.0 ) * ws ( k )
!Removed by Zhixiao. Water vapor mixing ratio (w) is not conserved during parcel lifting processes. We should not use w to define MU layer.
!thetav(k) = Theta ( VirtualTemperature ( tK (k), w (k) ), p(k)/100.0 )
!Added by Zhixiao. Critical modification: We use the model level with maximum equivalent potential temperature (etheta) below 500mb to define the MU layer
!Because equivalent potential temperature is conserved in dry and moist adiabatic processes.
etheta(k) = Thetae( tK(k), p(k)/100.0, rh(k), w(k) )
END DO
srclev = sfc
Expand All @@ -558,32 +570,38 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
srcws = ws ( sfc )
srcw = w ( sfc )
srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
srce = etheta (sfc) !Modified by Zhixiao
!~ 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
!Change initial searching level from the second to first model level. Because we did not compute pdiff, and p(k-1) properties is unnecessary.
!Modified by Zhixiao.
DO k = sfc, 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 we've made it past the lowest 500 hPa, exit the loop. MU layer is assumed below 500 hPa. Modified by Zhixiao.
! -------------------------------------------------------
IF ( pdiff >= REAL (180) ) EXIT
IF ( p(k) <= REAL (50000) ) EXIT

IF ( prcl == 1 ) THEN
! Removed by Zhixiao, w can not used for defining MU layer
!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
IF ( (w(k) > srcw) ) THEN
! Modified by Zhixiao, MU layer is featured by the max etheta
IF (etheta(k) > srce) THEN !Modified by Zhixiao.
srctheta = Theta ( tK(k), p(k)/100.0 )
srcw = w ( k )
srclev = k
srctK = tK ( k )
srcrh = rh ( k )
srcp = p ( k )
srce = etheta(k) !Modified by Zhixiao
END IF
END IF

Expand All @@ -605,12 +623,22 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
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) )


!Add caculation for LCL. Modified by Zhixiao
lcllev = -1 !Modified by Zhixiao
flag=.false. !Modified by Zhixiao
DO k=sfc,nz !zhixiao search the layer of LCL
IF (p (k) <= plcl) THEN !Modified by Zhixiao
lcllev=k !Modified by Zhixiao
flag=.true. !Modified by Zhixiao
END IF !Modified by Zhixiao
IF (flag) EXIT !Modified by Zhixiao
END DO !Modified by Zhixiao
flag=.false. !Modified by Zhixiao
!~ Now lift the parcel
! -------------------

Expand Down Expand Up @@ -728,47 +756,34 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
! -----------------------------------
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
ellev = -1 !Modified by Zhixiao
DO k = sfc, nz !Modified by Zhixiao
!~ LFC is defiend as the highest level when negative buyancy turns postive.
! Let us ignore the lower LFCs, and keep the highest LFC as the final result
! -----------------------------------
IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) <= plcl ) THEN !Modified by Zhixiao
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.
!~ Take the Highest EL as final result. Modified by Zhixiao
! ----------------------------------------------------------------
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
IF (k >= 2) THEN !Modified by Zhixiao
IF ( flag .and. buoy (k) < REAL (0) .and. buoy (k-1) >= REAL (0)) THEN !Modified by Zhixiao
ellev = k !Modified by Zhixiao
END IF
END IF
END DO
IF (ellev >= 0) THEN !Modified by Zhixiao
pel = p(ellev) !Modified by Zhixiao
CIN=REAL ( 0 ) !Modified by Zhixiao
DO k = sfc+1, nz !Modified by Zhixiao
! CAPE and CIN is defined as integrated positive and negative buoyant energy between LCL and EL, respectively. Modified by Zhixiao
IF ( p (k) <= plcl .and. p (k) >= pel) THEN !Modified by Zhixiao
CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) ) !Modified by Zhixiao
CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) ) !Modified by Zhixiao
END IF !Modified by Zhixiao
END DO !Modified by Zhixiao
END IF !Modified by Zhixiao
!~ Calculate lifted index by interpolating difference between
!~ parcel and ambient Tv to 500mb.
! ----------------------------------------------------------
Expand All @@ -783,7 +798,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
IF ( pd .le. pm ) THEN
lidx = 0.
EXIT

ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN

!~ Found trapping pressure: up, middle, down.
Expand All @@ -797,22 +812,21 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, &
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 )

IF ( CIN /= CIN ) cin = RUNDEF

!~ Chirp
! -----
Expand Down

0 comments on commit 9662d57

Please sign in to comment.