From 0d04d0d7c1d5a1549d3df89a00d2f92aad88ef48 Mon Sep 17 00:00:00 2001 From: hsoh-u Date: Fri, 18 Feb 2022 16:31:06 -0700 Subject: [PATCH] Feature 1824 pb2nc mlcape (#2057) Co-authored-by: Howard Soh --- met/data/config/PB2NCConfig_default | 3 +- met/docs/Users_Guide/reformat_point.rst | 3 +- met/src/tools/other/pb2nc/calcape.f | 28 ----- met/src/tools/other/pb2nc/pb2nc.cc | 147 ++++++++++++++++-------- test/config/PB2NCConfig_pbl | 3 +- 5 files changed, 105 insertions(+), 79 deletions(-) diff --git a/met/data/config/PB2NCConfig_default b/met/data/config/PB2NCConfig_default index e11fa94e30..13110d221d 100644 --- a/met/data/config/PB2NCConfig_default +++ b/met/data/config/PB2NCConfig_default @@ -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"; } ]; //////////////////////////////////////////////////////////////////////////////// diff --git a/met/docs/Users_Guide/reformat_point.rst b/met/docs/Users_Guide/reformat_point.rst index 7a3fce768f..0953afdd7c 100644 --- a/met/docs/Users_Guide/reformat_point.rst +++ b/met/docs/Users_Guide/reformat_point.rst @@ -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. _____________________ @@ -248,6 +248,7 @@ _____________________ { key = 'D_PRMSL'; val = 'PRMSL'; }, { key = 'D_PBL'; val = 'PBL'; }, { key = 'D_CAPE'; val = 'CAPE'; } + { key = 'D_MLCAPE'; val = 'MLCAPE'; } ]; diff --git a/met/src/tools/other/pb2nc/calcape.f b/met/src/tools/other/pb2nc/calcape.f index e8d3a60b9d..d833246671 100644 --- a/met/src/tools/other/pb2nc/calcape.f +++ b/met/src/tools/other/pb2nc/calcape.f @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 diff --git a/met/src/tools/other/pb2nc/pb2nc.cc b/met/src/tools/other/pb2nc/pb2nc.cc index cff540d107..0b914b3198 100644 --- a/met/src/tools/other/pb2nc/pb2nc.cc +++ b/met/src/tools/other/pb2nc/pb2nc.cc @@ -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]; @@ -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++) { @@ -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; @@ -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; @@ -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; @@ -1327,7 +1338,7 @@ void process_pbfile(int i_pb) { } } - if (cal_cape) { + if (cal_cape || cal_mlcape) { cape_level = 0; } @@ -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]; @@ -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++; } @@ -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; @@ -1646,14 +1657,6 @@ 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" ; @@ -1661,37 +1664,81 @@ void process_pbfile(int i_pb) { 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++; @@ -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"; } @@ -3003,7 +3050,11 @@ float compute_pbl(map 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 diff --git a/test/config/PB2NCConfig_pbl b/test/config/PB2NCConfig_pbl index 924c06a54c..ac48b5850f 100644 --- a/test/config/PB2NCConfig_pbl +++ b/test/config/PB2NCConfig_pbl @@ -92,7 +92,7 @@ level_category = [0, 1, 4, 5, 6]; // Use obs_bufr_map to rename variables in the output. // If empty, process all available variables. // -obs_bufr_var = ["D_CAPE", "D_PBL"]; +obs_bufr_var = ["D_CAPE", "D_PBL", "D_MLCAPE"]; //////////////////////////////////////////////////////////////////////////////// // @@ -120,6 +120,7 @@ obs_prepbufr_map = [ { key = "D_PBL"; val = "HPBL"; }, { key = "D_PRMSL"; val = "PRMSL"; }, { key = "D_CAPE"; val = "CAPE"; }, + { key = "D_MLCAPE"; val = "MLCAPE"; }, { key = "TDO"; val = "DPT"; }, { key = "PMO"; val = "PRMSL"; }, { key = "TOCC"; val = "TCDC"; },