Skip to content

Commit

Permalink
Feature 1824 pb2nc mlcape (#2057)
Browse files Browse the repository at this point in the history 
Co-authored-by: Howard Soh <[email protected]>
  • Loading branch information
@hsoh-u
hsoh-u and Howard Soh authored Feb 18, 2022
1 parent 418ebe7 commit 0d04d0d
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 79 deletions.
3 changes: 2 additions & 1 deletion met/data/config/PB2NCConfig_default
View file
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
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
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
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

0 comments on commit 0d04d0d

Please sign in to comment.