program pcpconform ! TAKEN FROM VERF_PRECIP FOR USE IN EVS GLOBAL_DET ATMOS ! ! the program reads in a precip grib file and convert it to 'regular' ! GRIB format, i.e. parameter is No. 61 (APCP) in GRIB1 Table 2. ! parameter(lmax=10000000) ! dimension f1(lmax), f2(lmax) integer ipopt(20), jpds(25), jgds(22), kpds(25), kgds(22) logical*1 bit(lmax) character arg*256,file1*200, file2*200, model*20 ! call getarg(1,arg) model=trim(arg) lmodel=len_trim(model) ! call getarg(2,arg) file1=trim(arg) len1=len_trim(file1) ! call getarg(3,arg) file2=trim(arg) len2=len_trim(file2) ! ! read in the input precip file using getgb: ! write(6,*) 'file1=', file1(1:len1) call baopenr(11,file1(1:len1),ibaret) call baopenw(51,file2,ibaret) ! jpds=-1 call getgb(11,0,lmax,-1,jpds,jgds,kf,k,kpds,kgds,bit,f1,iret) write(6,*) 'finished getgb, ibaret=', ibaret, ' iret=', iret ! if (model(1:3).eq.'gfs') then if (kpds(5).ne.59) then write(6,*) 'KPDS(5)=', kpds(5),' for ',model, ' STOP!' stop endif ! kpds(5)=61 dt = float(kpds(15)-kpds(14))*3600. f2 = f1*dt ! Change time range indicator from '3' (average) to '4' (accumulated ! from t1 to t2): kpds(16)=4 ! The original GFS precip rate (kg/m^2/s) has a Decimal Scaling factor of 6. ! This is way too high when we're converting the value to APCP. When testing ! the (3072 x 1536) GFS para output (Jan-Feb 2014), this leads to occasional ! errors in the output (no errors shown by putgb, but "missing end section ..." ! msg from wgrib and "GRIB header not found in scanning between records ..." ! from gribscan. kpds(22)=3 elseif (model(1:3).eq.'dwd') then ! 2014/7/1 ! DWD precip changed from GRIB1 to GRIB2 on 2014/6/25, but the GRIB2 parm ! is not 'precip'. Coping with it by converting it to GRIB1, then change ! it to the standard 'Table 2, parm 61' for precip. ! Also, after 'cnvgrib -g21', the DWD data somehow has decimal scale factor ! of 0 (even though the field - before I converted it to 'Table 2, parm 61' - ! seemed to have sufficient significant digits after decimal point: ! "min/max data 0 342.406"), and after resetting the grib1 file to ! 'Table 2, parm 61', it had no significant digits after the decimal point ! "min/max data 0 342"). So, set the decimal scaling factor to 3. ! ! 2015/1/26 DWD time info needs special treatment. Get the starting/ending ! hour from the last seven characters of the file name ! ! 2016/10/07 DWD GRIB2 changed again (see verflog.txt for details). Use ! cnvgrib2to1.py instead of cnvgrib for conversion. Still need some ! tweaking in pcpconform: ! 1) After cnvgrib2to1, kpds(13)=12: the forecast time unit is '12 hours'. ! this is inconsistent with kpds(14) and kpds(15): ! e.g., for dwd_2016100500_012_036.cnvgrib2to1, ! kpds(14)=60, kpds(15)=180, should have been '1' and '3' if the unit were ! indeed '12 hours'. Seems a bit dangerous simply devide kpds(14/15) by ! '60' - what if python library changes in the future in an unpredicted ! way? I'll use the time label in the file name and set time unit as ! 'hour', as in 2015/1/26. ! 2) Change longitude range from [0,-0.25] to [0,359.75] ! 3) Some small negative precip values - set min value to 0. ! 4) Still need to set the decimal scale factor to 3 - other wise the ! dec scaling factor is 0 and the field has no significant digits after ! the decimal point. ! ! 2016/10/25 pygrib not in production and too hard to have it installed in ! production. Revert to using cnvgrib. ! ! 1) Change parameter table to version 2, parameter to 61 (APCP) ! 1) Use the time label in the file name and set time unit as ! 'hour', as in 2015/1/26. ! 2) Change longitude range from [0,-0.25] to [0,359.75] ! 3) Some small negative precip values - set min value to 0. ! 4) Still need to set the decimal scale factor to 3 - other wise the ! dec scaling factor is 0 and the field has no significant digits after ! the decimal point. ! ! file1=dwd_2016100500_012_036.tmp ! ----:----|----:----|----:----| read(file1(len1-10:len1-4),'(i3,x,i3)') kpds(14), kpds(15) kpds(5)=61 ! parameter: APCP (was '255') kpds(13)=1 ! forecast time unit is 'hour' (was '0') kpds(16)=4 ! forecast time indiciator accumulation kpds(19)=2 ! Version number of parameter table (was '255') kpds(22)=3 ! Decimal scale factor (was '0') kgds(8)=359750 ! Longitude of the extreme point; was '-2147233) f2=max(0.,f1) elseif (model(1:5).eq.'ecmwf') then if (kpds(5).ne.228) then write(6,*) 'KPDS(5)=', kpds(5),' for ',model, ' STOP!' endif ! kpds(5)=61 ! Flip time range #1 and #2: kpds(15)=kpds(14) kpds(14)=0 ! Change time range indicator from '0' (instantaneous) to '4' (accumulated ! from t1 to t2): kpds(16)=4 ! Flip latitude range #1 and #2: kgdsdum=kgds(4) kgds(4)=kgds(7) kgds(7)=kgdsdum ! ! Flip the array in the north-south direction; also change precip ! accumulation unit from m to mm (i.e. kg/m2): nx=kgds(2) ny=kgds(3) k2 = 0 do 20 j = 1, ny do 10 i = 1, nx k2 = k2 + 1 k1 = (ny-j)*nx + i f2(k2) = f1(k1)*1000. 10 continue 20 continue ! ! Change PDS(1) from 98 (originating center: ECMWF) to 07 (NCEP). It seems ! that GEMPAK is hard-wired to flip the ECMWF plots, since the original ! ECMWF data are oriented north-to-south. kpds(1)=7 ! ! 'lspa' is #154 on Table 130. Convert to #61 on Table 2. elseif (model(lmodel-3:lmodel).eq.'soil') then kpds(5)=61 kpds(19)=2 f2=f1 endif ! call putgb(51,kf,kpds,kgds,bit,f2,iret) write(6,*) 'finished putgb, iret=', iret ! close(51) ! stop end