Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Feature 1824 pb2nc mlcape #2057

Merged
merged 10 commits into from
Feb 18, 2022
Merged

Feature 1824 pb2nc mlcape #2057

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
ef7b57f
#1824 Added MLCAPE
Feb 3, 2022
a6ba1ab
#1824 Cleanup
Feb 3, 2022
5dffde2
#1824 Support mixed layer cape
Feb 3, 2022
254ebf7
#1824 Support mixed layer cape
Feb 3, 2022
30f42d9
Merge remote-tracking branch 'origin/develop' into feature_1824_pb2nc…
Feb 17, 2022
432f7e0
#1824 ci-run-test Apply ivirt 1 and itype 2 only for MLCAPE and cleanup
Feb 18, 2022
656a103
#1824 Added MLCAPE to obs_bufr_map
Feb 18, 2022
bdf70b7
#1824 ci-run-test FIx for the Docker build
Feb 18, 2022
f8d6888
#1824 Added MLCAPE output into the existing unit test pb2nc_compute_p…
Feb 18, 2022
d44aa41
Merge remote-tracking branch 'origin/develop' into feature_1824_pb2nc…
Feb 18, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion met/data/config/PB2NCConfig_default
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,8 @@ obs_prepbufr_map = [
{ key = "D_MIXR"; val = "MIXR"; },
{ key = "D_PRMSL"; val = "PRMSL"; },
{ key = "D_PBL"; val = "PBL"; },
{ key = "D_CAPE"; val = "CAPE"; }
{ key = "D_CAPE"; val = "CAPE"; },
{ key = "D_MLCAPE"; val = "MLCAPE"; }
];

////////////////////////////////////////////////////////////////////////////////
Expand Down
3 changes: 2 additions & 1 deletion met/docs/Users_Guide/reformat_point.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ _____________________
obs_bufr_var = [ 'QOB', 'TOB', 'ZOB', 'UOB', 'VOB' ];


Each PrepBUFR message will likely contain multiple observation variables. The **obs_bufr_var** variable is used to specify which observation variables should be retained or derived. The variable name comes from BUFR file which includes BUFR table. The following BUFR names may be retained: QOB, TOB, ZOB, UOB, and VOB for specific humidity, temperature, height, and the u and v components of winds. The following BUFR names may be derived: D_DPT, D_WIND, D_RH, D_MIXR, D_PRMSL, D_PBL, and D_CAPE for dew point, wind speed, relative humidity, mixing ratio, pressure reduced to MSL, planetary boundary layer height, and convective available potential energy. This configuration replaces **obs_grib_code**. If the list is empty, all BUFR variables are retained.
Each PrepBUFR message will likely contain multiple observation variables. The **obs_bufr_var** variable is used to specify which observation variables should be retained or derived. The variable name comes from BUFR file which includes BUFR table. The following BUFR names may be retained: QOB, TOB, ZOB, UOB, and VOB for specific humidity, temperature, height, and the u and v components of winds. The following BUFR names may be derived: D_DPT, D_WIND, D_RH, D_MIXR, D_PRMSL, D_PBL, D_CAPE, and D_MLCAPE for dew point, wind speed, relative humidity, mixing ratio, pressure reduced to MSL, planetary boundary layer height, convective available potential energy, and mixed layer convective available potential energy. This configuration replaces **obs_grib_code**. If the list is empty, all BUFR variables are retained.

_____________________

Expand All @@ -248,6 +248,7 @@ _____________________
{ key = 'D_PRMSL'; val = 'PRMSL'; },
{ key = 'D_PBL'; val = 'PBL'; },
{ key = 'D_CAPE'; val = 'CAPE'; }
{ key = 'D_MLCAPE'; val = 'MLCAPE'; }
];


Expand Down
28 changes: 0 additions & 28 deletions met/src/tools/other/pb2nc/calcape.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,6 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
CINS(I,J) = 0.0
THESP(I,J)= 0.0
IEQL(I,J) = LM+1
C print*,' DEBUG HS calcape IEQL(I,J): by LM+1: ',IEQL(I,J)
C print*,' DEBUG HS calcape LCL(I,J) to 0 '

LCL(I,J)=0
ENDDO
Expand Down Expand Up @@ -233,11 +231,8 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
APESPK=(H10E5/TPSPK)**CAPA
TTHESK=TTHBTK*EXP(ELOCP*QBTK*APESPK/TTHBTK)
C--------------CHECK FOR MAXIMUM THETA E--------------------------------
C print*,' DEBUG HS calcape PTBL(IQ ,IT )',PTBL(IQ ,IT ),
C 2 PTBL(IQ+1 ,IT ),PTBL(IQ ,IT+1 ),PTBL(IQ+1 ,IT+1 )
C print*,' DEBUG HS calcape TTHESK.GT.THESP(I,J:',TTHESK,THESP(I,J)
IF(TTHESK.GT.THESP(I,J)) THEN
C print*,' DEBUG HS calcape Update PSP:',TPSPK,' LMM:',LMM
PSP (I,J)=TPSPK
THESP(I,J)=TTHESK
ENDIF
Expand All @@ -253,9 +248,7 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
DO J=1,JM
DO I=1,IM
IF (L.LT.LMH(I,J)) THEN
C print*,' DEBUG HS calcape L.LT.LMH(I,J):',L,LMH(I,J)
PKL = P(I,J,L)
C print*,' DEBUG HS calcape L,PKL.LT.PSP(I,J)',L,PKL,PSP(I,J)
IF (PKL.LT.PSP(I,J)) THEN
LCL(I,J)=L+1
PLCL(I,J)=P(I,J,L+1)
Expand All @@ -276,7 +269,6 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
DO J=1,JM
DO I=1,IM
PKL = P(I,J,L)
C print*,' DEBUG HS calcape L.LE.LCL(I,J):',L,LCL(I,J)
IF(L.LE.LCL(I,J)) THEN
IF(PKL.LT.PLQ)THEN
KNUML=KNUML+1
Expand Down Expand Up @@ -311,23 +303,19 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
ENDIF

C------------SEARCH FOR EQ LEVEL----------------------------------------
C print*,' DEBUG HS calcape SEARCH FOR EQ LEVEL, KNUMH:',KNUMH
DO N=1,KNUMH
I=IHRES(N)
J=JHRES(N)
IF(TPAR(I,J,L).GT.T(I,J,L)) THEN
IEQL(I,J)=L
C print*,' DEBUG HS calcape 111 IEQL(I,J):',IEQL(I,J)
PEQL(I,J)=P(I,J,L)
ENDIF
ENDDO
C print*,' DEBUG HS calcape SEARCH FOR EQ LEVEL, KNUML:',KNUML
DO N=1,KNUML
I=ILRES(N)
J=JLRES(N)
IF(TPAR(I,J,L).GT.T(I,J,L)) THEN
IEQL(I,J)=L
C print*,' DEBUG HS calcape 222 IEQL(I,J):',IEQL(I,J)
PEQL(I,J)=P(I,J,L)
ENDIF
ENDDO
Expand All @@ -339,7 +327,6 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
DO I=1,IM
LCLK=LCL(I,J)
IEQK=IEQL(I,J)
C print*,' DEBUG HS calcape IEQK,LCLK:',IEQK,LCLK
DO L=IEQK,LCLK
c print*,'l=',l
c print*,'p(i,j,l)=',p(i,j,l)
Expand Down Expand Up @@ -369,12 +356,10 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
C print*,'virtual thetap=',thetap
C print*,'virtual thetaa=',thetaa
endif
C print*,' DEBUG HS calcape before CAPE(I,J)=',CAPE(I,J)
IF (THETAP.LT.THETAA)
& CINS(I,J)=CINS(I,J)+G*(ALOG(THETAP)-ALOG(THETAA))*DZKL
IF (THETAP.GT.THETAA)
& CAPE(I,J)=CAPE(I,J)+G*(ALOG(THETAP)-ALOG(THETAA))*DZKL
C print*,' DEBUG HS calcape after CAPE(I,J)=',CAPE(I,J)
ENDDO
ENDDO
ENDDO
Expand All @@ -389,7 +374,6 @@ SUBROUTINE CALCAPE(ivirt,itype,T,Q,P,p1d,t1d,q1d,PINT,LMH,
DO I = 1,IM
CAPE(I,J) = AMAX1(0.0,CAPE(I,J))
CINS(I,J) = AMIN1(CINS(I,J),0.0)
C print*,' DEBUG HS calcape final CAPE(I,J)=',CAPE(I,J)
ENDDO
ENDDO
c if(itype.eq.2) print*,'CAPE,CINS=',CAPE,CINS
Expand Down Expand Up @@ -791,19 +775,7 @@ SUBROUTINE SPLINE(JTB,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
600 K1=K1+1
IF(K1.LE.NNEW) GO TO 300
C-----------------------------------------------------------------------
C print*,' DEBUG HS calcape NOLD,NNEW: ',NOLD,NNEW
CC print*,' DEBUG HS calcape XOLD: ',XOLD
CC print*,' DEBUG HS calcape YOLD: ',YOLD
C print*,' DEBUG HS calcape XOLD(1)->: XNEW(1)',XOLD(1),XNEW(1)
C print*,' DEBUG HS calcape XOLD(2)->: XNEW(2)',XOLD(2),XNEW(2)
CC print*,' DEBUG HS calcape XNEW: ',XNEW
CC print*,' DEBUG HS calcape YNEW: ',YNEW
C print*,' DEBUG HS calcape YOLD(1)->: YNEW(1)',YOLD(1),YNEW(1)
C print*,' DEBUG HS calcape YOLD(2)->: YNEW(2)',YOLD(2),YNEW(2)
C print*,' DEBUG HS calcape YOLD(3)->: YNEW(3)',YOLD(3),YNEW(3)
C print*,' DEBUG HS calcape YOLD(-2)->: YNEW(-2)',
C 2 YOLD(NOLD-1),YNEW(NNEW-1)
C print*,' DEBUG HS calcape YOLD(-1)->: YNEW(-1)',
C 2 YOLD(NOLD),YNEW(NNEW)
RETURN
END
Expand Down
147 changes: 99 additions & 48 deletions met/src/tools/other/pb2nc/pb2nc.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -260,8 +260,10 @@ static const char *airnow_aux_vars = "TPHR QCIND";
// Pick the latter one if exists multiuple variables
static const char *bufr_avail_sid_names = "SID SAID RPID";
static const char *bufr_avail_latlon_names = "XOB CLON CLONH YOB CLAT CLATH";
static const char *derived_mlcape = "D_MLCAPE";
static const char *derived_cape = "D_CAPE";
static const char *derived_pbl = "D_PBL";
static const float MLCAPE_INTERVAL = 3000.;

static double bufr_obs[mxr8lv][mxr8pm];
static double bufr_obs_extra[mxr8lv][mxr8pm];
Expand Down Expand Up @@ -511,6 +513,7 @@ void initialize() {
prepbufr_derive_vars.add("D_MIXR");
prepbufr_derive_vars.add("D_PRMSL");
prepbufr_derive_vars.add("D_CAPE");
prepbufr_derive_vars.add("D_MLCAPE");
prepbufr_derive_vars.add("D_PBL");

for (int idx=0; idx<(sizeof(hdr) / sizeof(hdr[0])); idx++) {
Expand Down Expand Up @@ -982,10 +985,13 @@ void process_pbfile(int i_pb) {
int itype = 1; // itype 1: where a parcel is lifted from the ground
// itype 2: Where the "best cape" in a number of parcels
int cape_code = -1;
int mlcape_code = -1;
float p1d,t1d,q1d;
int IMM, JMM;
int cape_level=0, cape_count=0, cape_cnt_too_big=0, cape_cnt_surface_msgs=0;
int cape_cnt_no_levels=0, cape_cnt_missing_values=0, cape_cnt_zero_values=0;
int cape_level, cape_count, cape_cnt_too_big, cape_cnt_surface_msgs;
int cape_cnt_no_levels, cape_cnt_missing_values, cape_cnt_zero_values;
int mlcape_count, mlcape_cnt_too_big;
int mlcape_cnt_missing_values, mlcape_cnt_zero_values;
float cape_p, cape_h;
float cape_qm = bad_data_float;

Expand All @@ -997,8 +1003,9 @@ void process_pbfile(int i_pb) {

bool has_pbl_data;
bool do_pbl = false;
bool cal_cape = bufr_obs_name_arr.has(derived_cape, cape_code, false);
bool cal_pbl = bufr_obs_name_arr.has(derived_pbl, pbl_code, false);
bool cal_cape = bufr_obs_name_arr.has(derived_cape, cape_code, false);
bool cal_mlcape = bufr_obs_name_arr.has(derived_mlcape, mlcape_code, false);

bool is_same_header;
unixtime prev_hdr_vld_ut = (unixtime) 0;
Expand All @@ -1010,6 +1017,10 @@ void process_pbfile(int i_pb) {

// Initialize
prev_hdr_lat = prev_hdr_lon = prev_hdr_elv = bad_data_double;
cape_level = cape_count = cape_cnt_too_big = 0;
cape_cnt_no_levels = cape_cnt_missing_values = cape_cnt_zero_values = 0;
mlcape_count = mlcape_cnt_too_big = 0;
mlcape_cnt_missing_values = mlcape_cnt_zero_values = 0;

if (cal_pbl) {
is_same_header = false;
Expand Down Expand Up @@ -1327,7 +1338,7 @@ void process_pbfile(int i_pb) {
}
}

if (cal_cape) {
if (cal_cape || cal_mlcape) {
cape_level = 0;
}

Expand Down Expand Up @@ -1471,7 +1482,7 @@ void process_pbfile(int i_pb) {
obs_arr[4] += c_to_k;
}

if (cal_cape) {
if (cal_cape || cal_mlcape) {
if(grib_code == pres_grib_code) {
is_cape_input = true;
if (cape_level < MAX_CAPE_LEVEL) cape_data_pres[cape_level] = obs_arr[4];
Expand Down Expand Up @@ -1567,7 +1578,7 @@ void process_pbfile(int i_pb) {
} // end for i
}

if (cal_cape) {
if (cal_cape || cal_mlcape) {
if (cape_member_cnt >= 3) cape_level++;
}

Expand Down Expand Up @@ -1603,7 +1614,7 @@ void process_pbfile(int i_pb) {
}
} // end for lv

if (cal_cape) {
if (cal_cape || cal_mlcape) {
if (1 < cape_level) {
bool reverse_levels;
float cape_val, cin_val, PLCL,PEQL;
Expand Down Expand Up @@ -1646,52 +1657,88 @@ void process_pbfile(int i_pb) {
}
}

//p1d = cape_p;
//t1d = cape_data_temp[cape_level-1];
//q1d = cape_data_spfh[cape_level-1];
calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres,
&p1d,&t1d,&q1d, static_dummy_201,
&cape_level, &IMM,&JMM, &cape_level,
&cape_val, &cin_val, &PLCL, &PEQL, static_dummy_200);

if(mlog.verbosity_level() >= 7) {
mlog << Debug(7) << method_name << " index,P,T,Q to compute CAPE from "
<< i_read << "-th message\n" ;
for (int idx=0; idx<cape_level; idx++) {
mlog << Debug(7) << method_name << " " << idx << ", " << cape_data_pres[idx]
<< ", " << cape_data_temp[idx] << ", " << cape_data_spfh[idx] << "\n";
}
mlog << Debug(7) << method_name
<< " calcape_(" << ivirt << "," << itype << ") cape_val: "
<< cape_val << " cape_level: " << cape_level
<< ", cin_val: " << cin_val
<< ", PLCL: " << PLCL << ", PEQL: " << PEQL
<< " lat: " << hdr_lat << ", lon: " << hdr_lon
<< " valid_time: " << unix_to_yyyymmdd_hhmmss(hdr_vld_ut)
<< " " << hdr_typ << " " << hdr_sid
<< "\n\n" ;
}

if (cape_val > MAX_CAPE_VALUE) {
cape_cnt_too_big++;
mlog << Debug(5) << method_name
<< " Ignored cape_value: " << cape_val << " cape_level: " << cape_level
<< ", cin_val: " << cin_val
<< ", PLCL: " << PLCL << ", PEQL: " << PEQL << "\n";
if (cal_cape) {
calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres,
&p1d,&t1d,&q1d, static_dummy_201,
&cape_level, &IMM,&JMM, &cape_level,
&cape_val, &cin_val, &PLCL, &PEQL, static_dummy_200);

if(mlog.verbosity_level() >= 7) {
mlog << Debug(7) << method_name
<< " calcape_(" << ivirt << "," << itype << ") cape_val: "
<< cape_val << " cape_level: " << cape_level
<< ", cin_val: " << cin_val
<< ", PLCL: " << PLCL << ", PEQL: " << PEQL
<< " lat: " << hdr_lat << ", lon: " << hdr_lon
<< " valid_time: " << unix_to_yyyymmdd_hhmmss(hdr_vld_ut)
<< " " << hdr_typ << " " << hdr_sid
<< "\n\n" ;
}

if (cape_val > MAX_CAPE_VALUE) {
cape_cnt_too_big++;
mlog << Debug(5) << method_name
<< " Ignored cape_value: " << cape_val << " cape_level: " << cape_level
<< ", cin_val: " << cin_val
<< ", PLCL: " << PLCL << ", PEQL: " << PEQL << "\n";
}
else if (cape_val >= 0) {
obs_arr[1] = cape_code;
obs_arr[2] = cape_p;
obs_arr[3] = cape_h;
obs_arr[4] = cape_val; // observation value
addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut,
hdr_lat, hdr_lon, hdr_elv, cape_qm,
OBS_BUFFER_SIZE);
cape_count++;
n_derived_obs++;
if (is_eq(cape_val, 0.)) cape_cnt_zero_values++;
}
else cape_cnt_missing_values++;
}
else if (cape_val >= 0) {
obs_arr[1] = cape_code;
obs_arr[2] = cape_p;
obs_arr[3] = cape_h;
obs_arr[4] = cape_val; // observation value
addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut,
hdr_lat, hdr_lon, hdr_elv, cape_qm,
OBS_BUFFER_SIZE);
cape_count++;
n_derived_obs++;
if (is_eq(cape_val, 0.)) cape_cnt_zero_values++;

if (cal_mlcape) {
float mlcape_val = bad_data_float;

ivirt = 1;
itype = 2;
// The second last seems to be better than the average of last two or three
p1d = cape_data_pres[cape_level-2];
t1d = cape_data_temp[cape_level-2];
q1d = cape_data_spfh[cape_level-2];
calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres,
&p1d,&t1d,&q1d, static_dummy_201,
&cape_level, &IMM,&JMM, &cape_level,
&mlcape_val, &cin_val, &PLCL, &PEQL, static_dummy_200);
if (mlcape_val > MAX_CAPE_VALUE) {
mlcape_cnt_too_big++;
mlog << Debug(5) << method_name
<< " Ignored ML_cape: " << mlcape_val << " cape_level: "
<< cape_level << "\n";
}
else if (mlcape_val >= 0) {
obs_arr[1] = mlcape_code;
obs_arr[2] = cape_p;
obs_arr[3] = cape_h;
obs_arr[4] = mlcape_val; // observation value
addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut,
hdr_lat, hdr_lon, hdr_elv, cape_qm,
OBS_BUFFER_SIZE);
mlcape_count++;
n_derived_obs++;
if (is_eq(mlcape_val, 0.)) mlcape_cnt_zero_values++;
}
else mlcape_cnt_missing_values++;
}
else cape_cnt_missing_values++;
}
else if (1 < buf_nlev) cape_cnt_no_levels++;
else cape_cnt_surface_msgs++;
Expand Down Expand Up @@ -1948,13 +1995,13 @@ void process_pbfile(int i_pb) {
<< "Total observations retained or derived\t= "
<< (n_file_obs + n_derived_obs) << "\n";

if (cal_cape) {
mlog << Debug(3) << "\nDerived CAPE = " << cape_count
<< "\tZero = " << cape_cnt_zero_values
if (cal_cape || cal_mlcape) {
mlog << Debug(3) << "\nDerived CAPE = " << (cape_count + mlcape_count)
<< "\tZero = " << (cape_cnt_zero_values + mlcape_cnt_zero_values)
<< "\n\tnot derived: No cape inputs = " << (cape_cnt_no_levels)
<< "\tNo vertical levels = " << (cape_cnt_surface_msgs)
<< "\n\tfiltered: " << cape_cnt_missing_values << ", "
<< cape_cnt_too_big
<< "\n\tfiltered: " << (cape_cnt_missing_values + mlcape_cnt_missing_values) << ", "
<< (cape_cnt_too_big + mlcape_cnt_too_big)
<< "\n";
}

Expand Down Expand Up @@ -3003,7 +3050,11 @@ float compute_pbl(map<float, float*> pqtzuv_map_tq,

if (pbl_level <= 0) {
mlog << Debug(4) << method_name
<< "Skip CALPBL because an empty list after combining TQZ and UV\n";
<< "Skip CALPBL because of an empty list after combining TQZ and UV\n";
}
else if (pbl_level == 1) {
mlog << Debug(4) << method_name
<< "Skip CALPBL because of only one available record after combining TQZ and UV\n";
}
else {
// Order all observations by pressure from bottom to top
Expand Down
Loading