From 1ac3a8e8fea004b2ca36c6f127dc977b9f280f95 Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:58:57 -0700 Subject: [PATCH 1/9] Modified module_diag_functions.F 1. Correct the most unstable (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 that is conserved during the dry/moist adiabatic lifting processes. 2. Correct precipitable water (Pwat) definition from column accumulated 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. Supplement equilibrium level (EL) calculation, which is the highest model level when buoyancy turns from positive to negative. Save corrected MU layer height, LCL, EL as output variables. 4. Compute and save 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. --- phys/module_diag_functions.F | 171 +++++++++++++++++------------------ 1 file changed, 83 insertions(+), 88 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 41049ce3cc..994f3bd2be 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -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 @@ -436,8 +439,8 @@ END FUNCTION Pwat !~ to higher levels. ~! !~ ~! !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! - FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & - parcel ) result (ostat) + FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & + parcel ) result (ostat) !Add declaration of lcl, el, and MU layer. Modified by Zhixiao. IMPLICIT NONE @@ -449,9 +452,10 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & 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 ) :: zlcl !~ LCL Height ( gpm ) Modified by Zhixiao. REAL, INTENT ( OUT ) :: zlfc !~ LFC Height ( gpm ) - REAL, INTENT ( OUT ) :: plfc !~ LFC Pressure ( Pa ) - REAL, INTENT ( OUT ) :: lidx !~ Lifted index + REAL, INTENT ( OUT ) :: zel !~ EL Height ( gpm ) Modified by Zhixiao. + REAL, INTENT ( OUT ) :: zmu !~ Most unstable parcel height ( gpm ) Modified by Zhixiao. INTEGER :: ostat !~ Function return status !~ Nonzero is bad. @@ -463,10 +467,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 !~ EL pressure ( Pa ). Modified by Zhixiao. REAL :: nbuoy !~ Negative buoyancy REAL :: pbuoy !~ Positive buoyancy @@ -481,6 +487,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 ! ---------------- @@ -493,13 +500,16 @@ 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 LOGICAL :: flag !~ Dummy flag LOGICAL :: wflag !~ Saturation flag REAL :: freeze !~ Water loading multiplier - REAL :: pdiff !~ Pressure difference between levs + REAL :: pdiff !~ Pressure difference between levs + REAL :: hdiff !~ Height difference between levs. Modified by Zhixiao. REAL :: pm, pu, pd !~ Middle, upper, lower pressures REAL :: lidxu !~ Lifted index at upper level REAL :: lidxd !~ Lifted index at lower level @@ -513,15 +523,17 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & REAL :: g !~ Acceleration due to gravity PARAMETER ( g = 9.80665 ) !~ m s^-2 REAL :: RUNDEF - PARAMETER ( RUNDEF = -9.999E30 ) + PARAMETER ( RUNDEF = -999 ) !Specialized for CACTI. CACTI-Simulation default filling value is -999. Modified by Zhixiao. Please reset RUNDEF = -9.999E30 for general use. - !~ Initialize variables + !~ Initialize variables zhixiao ! -------------------- 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 + ZLCL = RUNDEF ! Modified by Zhixiao + ZEL = RUNDEF ! Modified by Zhixiao + ZMU = RUNDEF ! Modified by Zhixiao !~ Look for submitted parcel definition !~ 1 = Most unstable @@ -549,6 +561,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 conseved both in dry and moist adiabatic processes. + etheta(k) = Thetae( tK(k), p(k)/100.0, rh(k), w(k) ) END DO srclev = sfc @@ -558,32 +575,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 @@ -605,12 +628,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 ! ------------------- @@ -728,91 +761,53 @@ 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 - - !~ 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 (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 + !~ Covert LCL and LFC from pressure level to height. Modified by Zhixiao ! --------------------------------------------------- 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 ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN + ! ZLFC = REAL ( -1 ) + ! END IF + + IF ( lcllev > 0 ) THEN !Modified by Zhixiao + ZLCL = hgt ( lcllev ) !Modified by Zhixiao + END IF !Modified by Zhixiao + + IF ( ellev > 0 ) THEN !Modified by Zhixiao + ZEL = hgt ( ellev ) !Modified by Zhixiao + END IF !Modified by Zhixiao + + ZMU = hgt ( srclev ) !Modified by Zhixiao !~ Chirp ! ----- From 20fc1d0c4c454fb3132187f62dfd00f6d1a9129b Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Wed, 10 Mar 2021 22:05:31 -0700 Subject: [PATCH 2/9] Bug fix of MU layer and PWAT --- phys/module_diag_functions.F | 1 - 1 file changed, 1 deletion(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 994f3bd2be..1c7a5df56c 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -509,7 +509,6 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & LOGICAL :: wflag !~ Saturation flag REAL :: freeze !~ Water loading multiplier REAL :: pdiff !~ Pressure difference between levs - REAL :: hdiff !~ Height difference between levs. Modified by Zhixiao. REAL :: pm, pu, pd !~ Middle, upper, lower pressures REAL :: lidxu !~ Lifted index at upper level REAL :: lidxd !~ Lifted index at lower level From 08c113bfc2e6ec5772dfe437fb6a41f33df0b0bc Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 10:52:26 -0700 Subject: [PATCH 3/9] Minimize modifications --- phys/module_diag_functions.F | 45 ++++++++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 1c7a5df56c..57559c6b3b 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -439,8 +439,8 @@ END FUNCTION Pwat !~ to higher levels. ~! !~ ~! !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! - FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & - parcel ) result (ostat) !Add declaration of lcl, el, and MU layer. Modified by Zhixiao. + FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & + parcel ) result (ostat) IMPLICIT NONE @@ -454,6 +454,8 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & REAL, INTENT ( OUT ) :: cin !~ CIN ( J kg^-1 ) REAL, INTENT ( OUT ) :: zlcl !~ LCL Height ( gpm ) Modified by Zhixiao. REAL, INTENT ( OUT ) :: zlfc !~ LFC Height ( gpm ) + REAL, INTENT ( OUT ) :: plfc !~ LFC Pressure ( Pa ) + REAL, INTENT ( OUT ) :: lidx !~ Lifted index REAL, INTENT ( OUT ) :: zel !~ EL Height ( gpm ) Modified by Zhixiao. REAL, INTENT ( OUT ) :: zmu !~ Most unstable parcel height ( gpm ) Modified by Zhixiao. INTEGER :: ostat !~ Function return status @@ -508,7 +510,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & LOGICAL :: flag !~ Dummy flag LOGICAL :: wflag !~ Saturation flag REAL :: freeze !~ Water loading multiplier - REAL :: pdiff !~ Pressure difference between levs + 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 @@ -522,9 +524,9 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & REAL :: g !~ Acceleration due to gravity PARAMETER ( g = 9.80665 ) !~ m s^-2 REAL :: RUNDEF - PARAMETER ( RUNDEF = -999 ) !Specialized for CACTI. CACTI-Simulation default filling value is -999. Modified by Zhixiao. Please reset RUNDEF = -9.999E30 for general use. - - !~ Initialize variables zhixiao + PARAMETER ( RUNDEF = -9.999E30 ) + + !~ Initialize variables ! -------------------- ostat = 0 CAPE = REAL ( 0 ) @@ -788,9 +790,38 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlcl, zlfc, zel, zmu, & END IF !Modified by Zhixiao END DO !Modified by Zhixiao END IF !Modified by Zhixiao - !~ Covert LCL and LFC from pressure level to height. Modified by Zhixiao + !~ 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 + !~ Save LFC pressure and height. Modified by Zhixiao ! --------------------------------------------------- IF ( lfclev > 0 ) THEN + PLFC = p ( lfclev ) ZLFC = hgt ( lfclev ) END IF From 5070041a723ec62216bc7be263f49da383532c1b Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 10:59:41 -0700 Subject: [PATCH 4/9] Adjust spaces --- phys/module_diag_functions.F | 50 ++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 57559c6b3b..4646fdba23 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -791,33 +791,33 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & END DO !Modified by Zhixiao END IF !Modified by Zhixiao !~ Calculate lifted index by interpolating difference between - !~ parcel and ambient Tv to 500mb. - ! ---------------------------------------------------------- - DO k = sfc + 1, nz + !~ 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 - 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 + ENDIF - END DO + END DO !~ Save LFC pressure and height. Modified by Zhixiao ! --------------------------------------------------- IF ( lfclev > 0 ) THEN From 4da9f32ec986fb1695206e32fe5e24d9e00ba8d5 Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 16:05:40 -0700 Subject: [PATCH 5/9] Clean up irrelevant calculations --- phys/module_diag_functions.F | 35 +++++++++-------------------------- 1 file changed, 9 insertions(+), 26 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 4646fdba23..7813d1454d 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -440,7 +440,7 @@ END FUNCTION Pwat !~ ~! !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & - parcel ) result (ostat) + parcel ) result (ostat) IMPLICIT NONE @@ -452,12 +452,9 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & 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 ) :: zlcl !~ LCL Height ( gpm ) Modified by Zhixiao. REAL, INTENT ( OUT ) :: zlfc !~ LFC Height ( gpm ) REAL, INTENT ( OUT ) :: plfc !~ LFC Pressure ( Pa ) REAL, INTENT ( OUT ) :: lidx !~ Lifted index - REAL, INTENT ( OUT ) :: zel !~ EL Height ( gpm ) Modified by Zhixiao. - REAL, INTENT ( OUT ) :: zmu !~ Most unstable parcel height ( gpm ) Modified by Zhixiao. INTEGER :: ostat !~ Function return status !~ Nonzero is bad. @@ -474,7 +471,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & REAL :: buoy ( nz ) !~ Buoyancy REAL :: tlclK !~ LCL temperature ( K ) REAL :: plcl !~ LCL pressure ( Pa ) - REAL :: pel !~ EL pressure ( Pa ). Modified by Zhixiao. + REAL :: pel !~ Equilibrium pressure ( Pa ). Modified by Zhixiao. REAL :: nbuoy !~ Negative buoyancy REAL :: pbuoy !~ Positive buoyancy @@ -525,16 +522,12 @@ 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 = 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 - ZLCL = RUNDEF ! Modified by Zhixiao - ZEL = RUNDEF ! Modified by Zhixiao - ZMU = RUNDEF ! Modified by Zhixiao + PLFC = RUNDEF !~ Look for submitted parcel definition !~ 1 = Most unstable @@ -565,7 +558,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & !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 conseved both in dry and moist adiabatic processes. + !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 @@ -820,24 +813,14 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & END DO !~ Save LFC pressure and height. Modified by Zhixiao ! --------------------------------------------------- - IF ( lfclev > 0 ) THEN - PLFC = p ( lfclev ) - ZLFC = hgt ( lfclev ) + IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN + PLFC = REAL ( -1 ) + ZLFC = REAL ( -1 ) END IF - ! IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN - ! ZLFC = REAL ( -1 ) - ! END IF + IF ( CAPE /= CAPE ) cape = REAL ( 0 ) - IF ( lcllev > 0 ) THEN !Modified by Zhixiao - ZLCL = hgt ( lcllev ) !Modified by Zhixiao - END IF !Modified by Zhixiao - - IF ( ellev > 0 ) THEN !Modified by Zhixiao - ZEL = hgt ( ellev ) !Modified by Zhixiao - END IF !Modified by Zhixiao - - ZMU = hgt ( srclev ) !Modified by Zhixiao + IF ( CIN /= CIN ) cin = RUNDEF !~ Chirp ! ----- From a0074ae8ab65d99224f2ab2e661c8fe25423e1cc Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 16:09:45 -0700 Subject: [PATCH 6/9] Clean up irrelevant notes --- phys/module_diag_functions.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 7813d1454d..e20941f377 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -811,11 +811,11 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & ENDIF END DO - !~ Save LFC pressure and height. Modified by Zhixiao + !~ Assuming the the LFC is at a pressure level for now ! --------------------------------------------------- IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN - PLFC = REAL ( -1 ) - ZLFC = REAL ( -1 ) + PLFC = REAL ( -1 ) + ZLFC = REAL ( -1 ) END IF IF ( CAPE /= CAPE ) cape = REAL ( 0 ) From e9873f1934869edba621b122c2cb407c240e255d Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 16:11:37 -0700 Subject: [PATCH 7/9] Update module_diag_functions.F --- phys/module_diag_functions.F | 1 + 1 file changed, 1 insertion(+) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index e20941f377..ec4d593d69 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -527,6 +527,7 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & ostat = 0 CAPE = 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 !~ Look for submitted parcel definition From 24b4c43e7428eee8d171a98b9934bfb1922809f3 Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 17:35:42 -0700 Subject: [PATCH 8/9] Update module_diag_functions.F --- phys/module_diag_functions.F | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index ec4d593d69..2e6f584003 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -814,6 +814,11 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & 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 ) From f7ccb734fc50bb229a6bc84b04f16d8099b182e5 Mon Sep 17 00:00:00 2001 From: zhixiaozhang <72099269+zhixiaozhang@users.noreply.github.com> Date: Fri, 12 Mar 2021 17:36:48 -0700 Subject: [PATCH 9/9] Update module_diag_functions.F --- phys/module_diag_functions.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phys/module_diag_functions.F b/phys/module_diag_functions.F index 2e6f584003..5ec57adfbc 100644 --- a/phys/module_diag_functions.F +++ b/phys/module_diag_functions.F @@ -815,8 +815,8 @@ FUNCTION Buoyancy ( nz, tk, rh, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, & !~ Assuming the the LFC is at a pressure level for now ! --------------------------------------------------- IF ( lfclev > 0 ) THEN - PLFC = p ( lfclev ) - ZLFC = hgt ( lfclev ) + PLFC = p ( lfclev ) + ZLFC = hgt ( lfclev ) END IF IF ( PLFC /= PLFC .OR. PLFC < REAL (0) ) THEN