Enhance PB2NC to derive Mixed-Layer CAPE (MLCAPE).

JohnHalleyGotway commented 3 years ago

Describe the Enhancement

On 6/8/21, NOAA/EMC requested that the PB2NC tool be enhanced to derive additional variations of CAPE. They are most interested in mixed-layer CAPE (representing the average characteristics of the boundary layer). Mixed-layer is the "flavor" of CAPE that was recommended to be verified in our models, via the metrics workshop.

While PB2NC does currently derive CAPE, it is surface-based CAPE and its done by calling this Fortran sub-routine: https://github.com/dtcenter/MET/blob/main_v10.0/met/src/tools/other/pb2nc/calcape.f That routine was taken from the NOAA/EMC VSDB verification code, but no such code exists for mixed-layer CAPE. However, it is defined in UPP.

This task is to enhance PB2NC to derive the mixed-layer variation of CAPE. Consider some of the implementation options:

We could keep the scope limited and lift/adapt the code for deriving MLCAPE from the UPP source code. Should we keep that code in Fortran or re-write in C++?
We could expand the scope considerably and investigate a more direct interface to the UPP derivation routines. But this would add a new dependency on UPP.

An example of a UPP module to compute CAPE can be found at:

Time Estimate

Estimate the amount of work required here. Issues should represent approximately 1 to 3 days of work.

Sub-Issues

Consider breaking the enhancement down into sub-issues.

[ ] Add a checkbox for each sub-issue here.

Relevant Deadlines

needs to be done by 7/31/21 if possible

Funding Source

2792541

Define the Metadata

Assignee

[x] Select engineer(s) or no engineer required: not sure yet
[x] Select scientist(s) or no scientist required: Perry S

Labels

[x] Select component(s)
[x] Select priority
[x] Select requestor(s)

Projects and Milestone

[x] Select Repository and/or Organization level Project(s) or add alert: NEED PROJECT ASSIGNMENT label
[x] Select Milestone as the next official version or Future Versions

Define Related Issue(s)

Consider the impact to the other METplus components.

[x] METplus, MET, METdatadb, METviewer, METexpress, METcalcpy, METplotpy There should be no impacts to other METplus repositories.

Enhancement Checklist

See the METplus Workflow for details.

[ ] Complete the issue definition above, including the Time Estimate and Funding Source.
[ ] Fork this repository or create a branch of develop. Branch name: feature_<Issue Number>_<Description>
[ ] Complete the development and test your changes.
[ ] Add/update log messages for easier debugging.
[ ] Add/update unit tests.
[ ] Add/update documentation.
[ ] Push local changes to GitHub.
[ ] Submit a pull request to merge into develop. Pull request: feature <Issue Number> <Description>
[ ] Define the pull request metadata, as permissions allow. Select: Reviewer(s) and Linked issues Select: Repository level development cycle Project for the next official release Select: Milestone as the next official version
[ ] Iterate until the reviewer(s) accept and merge your changes.
[ ] Delete your fork or branch.
[ ] Close this issue.

JohnHalleyGotway commented 3 years ago

Checking the UPP source code, I see subroutines named CALCAPE and CALCAPE2 in this file: https://github.com/NOAA-EMC/EMC_post/blob/7df88065167f8e5fcecf161a3d8eea7fc44b5cd1/sorc/ncep_post.fd/UPP_PHYSICS.f

Tracing the logic back, I see the following:

RQSTFLD.F identifies variable id 582 as the "MIXED LAYER CAPE":

1798       DATA IFILV(582),AVBL(582),IQ(582),IS(582),AVBLGRB2(582)      &
1799      &                      /1,'MIXED LAYER CAPE    ',157,116,     &
1800      &                       'CAPE ON spec_pres_above_grnd'/

MISCLN.f says that if MIXED LAYER CAPE is requested, then itype = 2:

3093          IF(IGET(582)>0)THEN
3094            FIELD1=.TRUE.
3095          ENDIF
3096          IF(IGET(583)>0)THEN
3097            FIELD2=.TRUE.
3098          ENDIF
3099 
3100          IF(FIELD1.OR.FIELD2)THEN
3101            ITYPE = 2

MISCLN.f calls CAPCAPE with ITYPE as an input argument and sets EGRID1 as an output variable name:

3121            CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1,           &
3122                         EGRID2,EGRID3,EGRID4,EGRID5)

MISCLN.f stores the values from EGRID1 into a variable named MLCAPE:

3131                      GRID1(I,J) = EGRID1(I,J)
3132                      IF (SUBMODELNAME == 'RTMA') MLCAPE(I,J)=GRID1(I,J)

MISCLN.f writes GRID1 data to variable id 582:
```
3139                 fld_info(cfld)%ifld=IAVBLFLD(IGET(582))
3140                 fld_info(cfld)%lvl=LVLSXML(1,IGET(582))
```
So presumably, MIXED LAYER CAPE is computed by calling the CALCAPE subroutine (from UPP_PHYSICS) with ITYPE=2. Those are the values we need to add here and store them to an observation named MLCAPE.

Note the following from RQSTFLD.F:

Surface based CAPE is variable ID 32.
BEST CAPE and CIN are variable ID's 566 and 567.
VTCAPE and VTCIN (for Virtual Temperature) are variable ID's 452 and 453.
MLCAPE and MLCIN (for Mixed Layer) are variable ID's 582 and 583 (call CALCAPE with itype = 2).
MUCAPE and MUCIN (for Most Unstable) are variable ID's 584 and 585 (call CALCAPE with itype = 1).
CAPE 0-3KM and CIN 0-3km for HRRR are variable ID's 950 and 951.
Downdraft CAPE for HRRR is variable ID 954.

This issue adds support for MLCAPE. @PerryShafran-NOAA what other variations should also be added?

JohnHalleyGotway commented 3 years ago

From Geoff M: My take is that all of the flavors of CAPE can potentially be addressed in a single update. We right now have 4 possible CAPE values in the UPP: surface-based, mixed-layer, most unstable lowest 180 mb, and most unstable lowest 300 mb. Right now, itype is set to 1 (for SBCAPE) which tells the code to search the raob report for the highest value of theta-e in the lowest 10 mb (set with DPBND). For the most unstable lowest 300 mb, we could set itype to 4 and assign DPBND=300.E2 and we'd be all set (after changing all itype.eq.1 references to itype.eq.1 .or. itype.eq.4) For the others, itype=2 could be the most unstable lowest 180 mb, and itype=3 could be the mixed layer (and change all itype.eq.2 references to itype.eq.2 .or. itype.eq.3). For these instances, the logic is already build into the code to deal with direct assignments of the parcel to lift (as opposed to searching a depth); it just isn't currently used. We would just have to write some code to compute the mixed layer and 180 mb most unstable parcels (they're computed in the post processor in a bit of an odd manner, but they come out as P1D, T1D, and Q1D) from the raob data.

hsoh-u commented 3 years ago

A message at BUFR contains 1D P,T,Q. How to make 2D P,T,Q or 3D PBND,TBND,QBND? Can anybody explain me im, jm, jsta and jend?

JohnHalleyGotway commented 3 years ago

@hsoh-u my assumption, or at least hope, was that we'd be able to use the EXISTING calcape.f code in PB2NC to derive these: https://github.com/dtcenter/MET/blob/main_v10.0/met/src/tools/other/pb2nc/calcape.f

We're calling calcape_ with the first 2 args = 1. I think we can just modify those args with different values to derive at least some of the variables. But I think that's the real task here, figuring out those details.

hertneky commented 3 years ago

@JohnHalleyGotway @TaraJensen @hsoh-u I was asked to provide UPP output for MLCAPE for 0-100 mb above ground. UPP only computes for the lowest 90 mb. I will note that it is quite odd with how it does this. While the user-defined config file does have 'level' parameters, it is a bit misleading as it doesn't matter what values you give it since it doesn't use the levels given.

CALCAPE uses an average of the 1st 3 layers of P, T, Q boundary layers (computed in BNDLYR.f) as the input. There are 6 layers total (0-30 mb, 30-60 mb, 60-90 mb, etc). By averaging the 1st 3, you get the 0-90 mb layer for MLCAPE.

I do have grib2 output for MLCAPE on cheyenne:/glade/work/hertneky/upp/upprt/postprd. This was a sever wx case from the UFS SRW and the grib2 files from UPP are NATLEV.036 and PRSLEV.036. The MLCAPE is defined by the level info in the files (90-0 mb above ground).

hertneky commented 3 years ago

@JohnHalleyGotway @hsoh-u My apologies, I had changed my testing directory name to /glade/work/hertneky/upp/testing/postprd, where the files can be found.

hsoh-u commented 3 years ago

There is only one MLCAPE value at GFSPRS.006:

913:2532956579:d=2020020400:CAPE:90-0 mb above ground:6 hour fcst:

GFSPRS.006 includes P, T, and Q input, right?

hertneky commented 3 years ago

@hsoh-u I had originally directed you to the LAM output (NATLEV.036/PRSLEV.036) since they are smaller in size, both for input and output, but the GFS output is fine as well if you don't have issue with the large file sizes.

Model input for UPP for the GFS run can be found in /glade/work/hertneky/upp/testing/data/gfs, which will include the T, P, Q vars that UPP ingests to process MLCAPE (specifically atmf006.nc). q -> spfh t -> tmp p -> pfull (or UPP computes 3D pressure from pressfc+dpres and then averages from 2 layers and uses that throughout)

While the current grib2 UPP output I provided does have these variables, they are interpolated to pressure levels, compared to the raw model input. If it is better, I can output T, SPFH, PRES on the model levels in the UPP output, instead of using it from the model data. Please let me know, it isn't hard to add.

You are correct, there is only one MLCAPE variable in the UPP output and is for 0-90 mb as mentioned previously.

hsoh-u commented 3 years ago

T, SPFH, PRES from you is better than collecting them by me to reduce a humane error. MET does not interpolate to pressure levels. Please collect T, SPFH, PRES which are interpolated to pressure levels) and give the P, T Q values to me.

hertneky commented 3 years ago

So the /glade/work/hertneky/upp/upprt/postprd/GFSPRS.006.dev has TMP, SPFH on isobaric levels, but I would think the data on model levels would be what you wanted since that's what UPP uses to compute MLCAPE?

hertneky commented 3 years ago

Perhaps I am just confused - I was under the impression you were trying to mimic the UPP algorithm for computing MLCAPE.

hsoh-u commented 3 years ago

Checked three cases but they did not match (MetPy failed to compute ML_CAPE with isobaric pressure levels)

ML-CAPE from UPP = 1451:

cape_val (1,1) = cape_val (1,2) = 1610.25
cape_val (0,1) = cape_val (0,2) = 1373.33

ML-CAPE from UPP = 1467:

cape_val(1,1) = cape_val(1,2) = 1607.9
cape_val(0,1) = cape_val(0,2) = 1369.98

ML-CAPE from UPP = 2232:

cape_val(1,1) = cape_val(1,2) = 2326.23
cape_val(0,1) = cape_val(0,2) = 2053.18

hertneky commented 3 years ago

@hsoh-u I would not expect the values to compare identically to UPP for these reasons: 1) UPPs computation is a bit convoluted in my opinion - see my previous mention on this. 2) UPPs MLCAPE is only for 90-0 mb above ground (not sure what you used for your tests) 3) UPP uses the values on model surfaces as input into the computation, not on isobaric surfaces - not sure how much this would change things, but imagine it would some.

JohnHalleyGotway commented 3 years ago

Met with @fossell and @hertneky on 10/7/21 to discuss.

There are 3 elements to this:

The pressure levels of the observations to be used as input.
The settings for ivirt (0 for regular theta, 1 for virtual theta) and itype (1 for surface-based parcel, 2 for best cape in the number of parcels)
What name to assign the result.

We need a table that consists of 4 columns: CAPE/CIN variation name, Input pressure levels, ivirt, itype

Recommend supporting all CAPE/CIN variations that are produced by UPP so that its output can be verified by observations. Ideally, we create a separate table for both UPP and calcape.f in MET to see how they compare.

See calcape.f: https://github.com/dtcenter/MET/blob/main_v10.0/met/src/tools/other/pb2nc/calcape.f

Plan:

@hertneky and @fossell document the CAPE/CIN UPP permutations in a table
@hsoh-u uses the table to make a best guess of those permutations in PB2NC
Commit the result, push feature branch to GitHub, compile on project machine(s)
Meet to plan logic for validation testing

hertneky commented 3 years ago

@JohnHalleyGotway @hsoh-u

I've created a document on the details of computing 4 different types of CAPE within UPP. Please note that the CAPE algorithm uses virtual temperature for all types (there is no IVIRT flag). Here is a link to the document:

https://docs.google.com/document/d/1IguVjWBX5vaQvZGfwfxrHqxovZoPHdZF1hNg3NFkDIs/edit

hsoh-u commented 2 years ago

from @hertneky:

 ITYPE=           2 DPBND=  0.0000000E+00 P1D(95,40)=   91592.19     T1D(95,40)=
   297.6417     Q1D(95,40)=  1.1071295E-02 MLCAPE(95,40)=   1282.799
 PMID(95,40,:)=   42.12350       101.0185       179.8740       270.1120
   373.3500       491.4290       626.4405       780.7574       957.0690
   1158.416       1388.231       1650.381       1949.208       2289.577
   2676.913       3117.249       3617.256       4184.272       4826.319
   5552.098       6370.962       7292.788       8327.331       9483.886
   10771.64       12199.70       13776.75       15510.72       17408.30
   19474.50       21712.09       24121.09       26698.39       29437.62
   32328.68       35357.30       38505.07       41749.65       45065.40
   48424.07       51795.77       55150.04       58456.98       61688.30
   64818.29       67824.55       70688.58       73396.05       75936.79
   78304.72       80497.49       82516.05       84364.09       86047.61
   87574.24       88952.89       90193.21       91305.30       92299.34
   93185.34       93973.05       94671.69       95290.00       95836.16
 T(95,40,:)=   258.8914       261.3713       260.2034       253.0265
   250.0544       245.9849       241.1160       236.4792       233.5308
   229.5910       225.9218       223.0825       222.0928       221.9948
   221.1809       218.5341       216.5348       215.7559       213.5920
   212.1257       211.9846       212.1231       210.9627       210.4264
   210.7383       211.8978       213.2945       214.2529       215.1741
   216.3933       219.1867       223.5339       227.8577       232.3464
   237.1857       242.0793       246.9203       251.6025       256.2554
   260.5364       263.8327       266.7532       270.9909       275.2119
   279.1491       282.7478       285.8876       288.4891       290.6883
   292.5164       294.4535       296.2153       297.4370       297.7391
   297.1093       296.6078       296.7967       297.3599       297.9618
   298.4389       298.8013       298.9513       298.9005       298.5872
 Q(95,40,:)=  2.5454005E-06  2.5354166E-06  2.4402493E-06  2.3346938E-06
  2.2826696E-06  2.2349423E-06  2.1807755E-06  2.1948333E-06  2.3504244E-06
  2.9591988E-06  4.1692465E-06  3.8182106E-06  3.4126410E-06  3.2589483E-06
  3.3125675E-06  3.4926888E-06  3.4944417E-06  2.9333878E-06  2.4024689E-06
  2.2679792E-06  2.4286846E-06  3.5141811E-06  4.2759366E-06  4.7172612E-06
  4.7531234E-06  5.6517019E-06  8.1814524E-06  1.2052872E-05  1.6188578E-05
  2.7972834E-05  3.2173346E-05  3.8350485E-05  6.2316984E-05  9.0379785E-05
  1.7427631E-04  2.5243388E-04  3.4471368E-04  3.9697401E-04  3.8653001E-04
  4.1821908E-04  9.1029552E-04  1.9419541E-03  2.3872368E-03  2.4001400E-03
  2.4071822E-03  2.5756054E-03  2.8330395E-03  3.5945713E-03  4.9012601E-03
  6.8181143E-03  7.8969793E-03  8.4347110E-03  8.5136052E-03  8.2985843E-03
  8.6999089E-03  9.7369943E-03  1.0827282E-02  1.1354864E-02  1.1779765E-02
  1.2168137E-02  1.2489138E-02  1.2814076E-02  1.3131263E-02  1.3479477E-02

Result by MET:

Was not able to compute the same p1d, t1d and q1d values.
By passing above p1d, t1d and q1d values:

Expected MLCAPE(95,40)=   1282.799
- itype should be 2

By using p1d, t1d, and q1d from the UPP

cape_val  (ivirt, itype)
   (1,1)=1936.41
   (1,2)=1295.92
   (0,1)=1829.93
   (0,2)=1212.57
   ( (0,2)+(1,2) ) / 2 =1254.24
   ( (0,1)+(1,1) ) / 2 =1883.17
   ( (0,1)+(0,2) ) / 2 =1521.25
   ( (1,1)+(1,2) ) / 2 =1616.16
   ( (1,1)+(0,2) ) / 2 =1574.49
   average of 4 values=1568.7

My MET (using second last)

Computed  p1d=91592.2 t1d=297.642 q1d=0.0110713
   (1,1)=1936.39
   (1,2)=1295.92
   (0,1)=1829.91
   (0,2)=1212.57
   (02&12)=1254.24
   ( (0,1) +(1,1) ) / 2=1883.15
   ( (0,1) +(0,2) ) / 2=1521.24  <-- MET's MLCAPE by using UPP's p1d, t1d, and q1d
   ( (1,1) +(1,2) ) / 2=1616.16
   ( (1,1) +(0,2) ) / 2=1574.48
   average of 4 values=1568.7

The computed MET's MLCAPE value is 1867.16 with p1d=95290 t1d=298.901 q1d=0.0131313 If a selected p1d is decreased, MLCAPE decreases, too.

Other cases (computing p1d instead of last second)

p1d=95563.1 t1d=298.744 q1d=0.0133054  AVG last two
cape_val (02&12)=1875.28

p1d=95265.9 t1d=298.813 q1d=0.0131416  AVG last three
cape_val (02&12)=1859.44

p1d=95836.2 t1d=298.587 q1d=0.0134795  last one
cape_val (02&12)=1883.15

hsoh-u commented 2 years ago

How to get above numbers:

Download the attached file (test_mlcape.txt) at pb2nc directory

cp -p pb2nc.cc pb2nc.cc.org
cat test_mlcape.txt >> pb2nc.cc

Open pb2nc.cc and rename the first main to main_org

 make
./pb2nc

test_mlcape.txt

hertneky commented 2 years ago

@hsoh-u

Result by MET:

* Was not able to compute the same p1d, t1d and q1d values.
* By passing above p1d, t1d and q1d values:
I'm not expecting values to be exact. From what I remember, there were some differences in the CALCAPE routine in UPP compared to what is used in MET. The values are close. I gave you input/output for MLCAPE from UPP which uses virtual temperature (equivalent to ivirt=1, itype=2 in MET). The result you got (1,2)=1295.92 is similar to UPP of 1282.779. If I am misunderstanding something here, then please let me know.
The computed MET's MLCAPE value is 1867.16 with p1d=95290 t1d=298.901 q1d=0.0131313
If a selected p1d is decreased, MLCAPE decreases, too.
I am fairly certain that the algorithm UPP uses to compute the T/P/Q of the parcel to be lifted is much different than what is in MET, in fact I don't know what MET does. UPP CALCAPE piggy-backs off another routine BNDLYR to get those means, which are mass weighted according to the documentation. Basically, this routine takes DPBND (specified layer) and computes the mean for a number of variables for that layer. Output from this routine is saved and then used in CALCAPE for MLCAPE.

hsoh-u commented 2 years ago

The MET's result were different for each cases (one method was better for one case, but not the other cases). If I have more sample data with expected MLCAPE values, the MET's logic will be improved.

hertneky commented 2 years ago

I'm not sure I understand, the 4 different computations you are providing are for the various types/methods for calculating CAPE. Is that what you meant by cases?

The MLCAPE inputs/outputs I gave is just one type of cape (mixed layer cape, hence why the parcel to be lifted is averaged over the 90 mb layer), so it is good that it is close in value to the MET MLCAPE value (itype=2, ivirt=1).

I sent you another data text file that also includes surface-based cape (itype=1, ivirt=1), for the same location, that you can compare to MET sfc cape. (It looks like MET sfc cape=1936.4.1 vs. 1913.8 from UPP, so again different but not too far fetched. In this data file - I also provided sfc cape and mixed layer cape at 2 other locations.

Note: I cannot provide data for comparing sfc cape and mixed layer cape for ivirt=0, since upp only uses virtual temperature. (corresponds to MET output for ivirt=0, itype=1 and ivirt=0, itype=2)

hsoh-u commented 2 years ago

Yes (x, y ) = (ivirt, itype)

Two more cases from @hertneky:

Case 2:

MET's MLCAPE
expected=1655.08  Using p1d=96953.6 t1d=296.207 q1d=0.0152234 (MET computed)
   (1,1)=1680.31
   (1,2)=1680.31
   (0,1)=1604.06
   (0,2)=1604.06
   (average of 02&12)=1642.18

expected=1655.08  Using p1d=91592.2 t1d=297.642 q1d=0.0110713 from UPP
   (1,1)=1680.31
   (1,2)=830.175
   (0,1)=1604.06
   (0,2)=784.365

Case 3:

MET's MLCAPE
expected=1755.88  Using  p1d=96514.7 t1d=294.806 q1d=0.0148598 (MET computed)
   (1,1)=1680.31
   (1,2)=1442.22
   (0,1)=1367.25
   (0,2)=1367.25
  (average of 02&12)=1664.7

expected=1755.88  Using pp1d=91592.2 t1d=297.642 q1d=0.0110713 from UPP
   (1,1)=1442.22
   (1,2)=985.928
   (0,1)=1367.25
   (0,2)=930.603

hsoh-u commented 2 years ago

Two more cases (the first one is the same as MLCAPE_output.txt cape_output.txt

TaraJensen commented 2 years ago

@hsoh-u - after reading through everything, I think it's okay that the values are not exactly what UPP computes due to what Tracy has stated several times. Based on the definition of mixed-layer CAPE from NOAA, the correct settings are ivrt=1 (virtual theta aka theta-e) and itype=2. Please go with those and push this into Pull Request with me as the reviewer.

hsoh-u commented 2 years ago

The unit test for PBL & CAPE were deleted. Is it OK to add an unit test for MLCAPE (bring back CAPE unit test)? Never mind. The test is moved to unit_pb2nc_indy.xml

JohnHalleyGotway commented 2 years ago

@hsoh-u I'm reopening this issue for a small amount of additional work requested.

During the NOAA METplus telecon on 3/7/2022, @PerryShafran-NOAA reported that his testing revealed that the derived MLCAPE values are reasonable, which is great news. However, he also noted that the log messages were a little confusing in that they failed to distinguish between CAPE and MLCAPE.

@TaraJensen requested that you create a feature branch to update the PB2NC log messages to differentiate between the derivation of CAPE and MLCAPE.

dtcenter / MET