------------------------------------------------------------------ This page is intended as a help file for using NCEP's Fortran GRIB decoding software to decode NCEP ensemble GRIB files. If you have any questions regarding anything on this page, please feel free to contact: timothy.marchok@noaa.gov yuejian.zhu@noaa.gov ------------------------------------------------------------------ Before getting started, there are 2 GRIB help files that I have online that you may find useful. If you are not familiar with the GRIB format, you should *definitely* read these help files before going on to the rest of this current help file, otherwise some of the material on the current page won't make much sense to you. http://www.emc.ncep.noaa.gov/gmb/tpm/grib/gribinfo.txt -- This is a basic intro describing what the GRIB format is. Read this one first. http://www.emc.ncep.noaa.gov/gmb/tpm/grib/gribinfo_for_w3lib.txt -- This I wrote for someone using our internal NCEP w3lib GRIB software, but if you're using the GRIB-reading software from NCEP, then this help file will apply to you as well, since the packages are very similar. The only caveat here is that I wrote this help file for someone who was reading regular GRIB files, NOT GRIB files that have an ensemble PDS extension. Therefore, the example that I show in that help file uses a call to "getgb" to read the data, while you will need a call to "getgbe" to read the data. The calls are very similar, but the getgbe call has 2 extra calling parameters (jens and kens) that you need in order to get the ensemble info that you need. Using the program included below, you can read GRIB records for whatever ensemble members you want. Before going on, definitely read the help file listed above, "gribinfo_for_w3lib.txt", specifically, the section titled "GETGB -- The /nwprod, GRIB-reading subroutine". Pay close attention to how you request what parameters you want to pull out of a file by specifying the array members jpds and jgds. Well, once you know how to do that, pulling out the particular ensemble member is done via the same method. You specify what member you want to read by setting both jens(2) and jens(3) to the needed values. Specifically, those values are as follows: For the 00z files: PDS Byte PDS Byte 42 43 jens(2) jens(3) Perturbation name (ens type) (ens ID) ----------------------- -------- ------ 1: High resolution control 1 1 2: Low resolution control 1 2 3: n1 2 1 4: p1 3 1 5: n2 2 2 6: p2 3 2 7: n3 2 3 8: p3 3 3 9: n4 2 4 10: p4 3 4 11: n5 2 5 12: p5 3 5 For the 12z files: PDS Byte PDS Byte 42 43 jens(2) jens(3) Perturbation name (ens type) (ens ID) ----------------------- -------- ------ 1: High resolution control 1 1 2: n1 2 1 3: p1 3 1 4: n2 2 2 5: p2 3 2 Note how in the program below, in the main loop, before each call to getgbe, I have the following code: jpds(23)=2 ! indicates that record has a PDS extension jens(2)=ens2ix(nf) jens(3)=ens3ix(nf) In this way, subroutine getgbe knows which ensemble member to search for in the GRIB file. By the way, note how in the program I have the following code: > if (nf.le.12) then > jpds(10)=gdd1 > jpds(11)=ghh1 > else > jpds(10)=gdd2 > jpds(11)=ghh2 > endif > jpds(14)=kftime(nt) I have this in here because I often look at ensemble files that have not only the current day's 00z data, but also the previous day's 12z data as well (you can "cat" GRIB files together on unix systems). From the ensemble member identifier table above, you can see that the 12z members will have the same jens(2) and jens(3) as some of the 00z members, so I have to request the specific forecast initial day and hour, as well as the forecast valid hour, in order to make sure that I'm getting the correct member. If you are reading files that have data from only the 00z run or only the 12z run, then you can eliminate the setting of those jpds(10) and jpds(11) parameters, and just leave the setting of the jpds(14) parameter so that you can get the forecast hour you want. In the program, I've also included 3 subroutines that you don't really need, but you can use them if you want. They are: subroutine conv1d2d -- getgbe returns a 1-d data array to you, and conv1d2d converts the data into a 2-d array. subroutine gribout -- Writes a GRIB record out to a new GRIB file. subroutine prtpds -- Prints the relevant pds and gds information out to standard output. Well, this should give you enough to get going. Let me know if you have any questions about any of this stuff or if you encounter any problems. Tim ______________________________________________________________________ Tim Marchok NCEP/Environmental Modeling Center Email: timothy.marchok@noaa.gov W/NP2, WWB Room 207 Phone: (301) 763-8000 x7264 Washington, DC 20233 Fax: (301) 763-8545 ______________________________________________________________________ --------------------------------------------------------------------------------- program cut_ensemble c c Example for unpacking ensemble grib data c c There is a 5 byte extension, starting at byte 41 in the pds, that identifies c the different ensemble members in grib data at NCEP. For other ensemble c products, such as precip probabilities, the extension goes beyond byte 45. c For more details on the extension, see: c c http://sgi62.wwb.noaa.gov:8080/ens/info/ens_grib.html#gribex c c EXTENSION OF PDS SECTION FOR THE IDENTIFICATION OF ENSEMBLE PRODUCTS c c OCTET IDENTIFIER CONTENT c ----- ---------- ------- c 41 Identifies application c 1 Ensemble c c 42 Type c 1 Unperturbed control forecast c 2 Individual negatively perturbed forecast c 3 Individual positively perturbed forecast c 4 Cluster c 5 Whole ensemble c c 43 Identification number c 1 (If byte 42=1, identifies high resolution control) c 2 (If byte 42=1, identifies low resolution control) c 3 c etc. c 44 Products c 1 Full field (individual forecast) / Unweighted mean c (cluster/ens) c 2 Weighted mean c 11 Standard deviation with respect to ensemble mean c 12 Standard deviation with respect to ensemble mean, c normalized c c 45 Spatial smoothing of product (number of highest total wavenumber c included) c 255 Original resolution retained c c Other change in main pds section: c c Octet 26 c c VALUE Interpretation c 2 NCEP Ensemble Products (A value of 2 in this byte (octet) c indicates that the current GRIB record contains an c extension to the PDS. That extension will contain c information relating to the identification of the c particular ensemble product. (See octets 41 and beyond) c By the way, the value in octet 26 will be returned to c your program from getgb or getgbe in kpds(23). c parameter(lugi=31,lugb=11,lugout=51,jf=512*256,nlat=73,nlon=144) character convflag*1,prtpdsflag*1,newgrbflag*1 integer jpds(25),jgds(22),kpds(25),kgds(22) integer jens(5),kens(5) integer ens2ix(17),ens3ix(17),kftime(33) logical lb(jf) real f(jf),data2d(nlon,nlat) namelist/namin/gdd1,ghh1,gdd2,ghh2 c data ens2ix/1,1,2,3,2,3,2,3,2,3,2,3,1,2,3,2,3/ data ens3ix/1,2,1,1,2,2,3,3,4,4,5,5,1,1,1,2,2/ data convflag/'y'/,prtpdsflag/'y'/,newgrbflag/'y'/ c read (5,namin,end=1000) 1000 continue c print *,'at beginning of tgetgbens.f, gdd1= ',gdd1,' ghh1= ',ghh1 & ,' gdd2= ',gdd2,' ghh2= ',ghh2 c j=0 c numfiles=17 numtimes=33 c do ik=1,numtimes kftime(ik) = (ik-1)*12 enddo c do nf=1,numfiles do nt=1,numtimes c print *,'--------------------------------------------' print *,'At top of loop, nf= ',nf,' nt= ',nt,' j=',j c Setting jpds and jgds tells subroutine getgb exactly what c record you want to read from the grib file. -1 is a wild card, c so that getgb would return the first record that it finds. So, c it's always a good idea to initially set both the jpds and jgds c arrays to -1, then set only the particular members of those c arrays in order for getgb to have enough information to find c your record. By the way, while jpds and jgds are the inputs to c getgb, kpds and kgds contain the output values, the actual c values, that getgb returns from its read of your grib file. c Similarly, set jens parameters to get the ensemble member you c want, and then getgbe returns the actual values that it found c in the record in kens. jpds=-1 jgds=-1 jens=-1 if (nf.le.12) then jpds(10)=gdd1 jpds(11)=ghh1 else jpds(10)=gdd2 jpds(11)=ghh2 endif jpds(14)=kftime(nt) c jpds(23)=2 ! indicates that record has a PDS extension jens(2)=ens2ix(nf) jens(3)=ens3ix(nf) c print *,'bgc, j= ',j,' nf= ',nf,' nt= ',nt & ,' jpds(2)= ',jpds(2),' jpds(3)= ',jpds(3) c call getgbe(lugb,lugi,jf,j,jpds,jgds,jens, & kf,k,kpds,kgds,kens,lb,f,iret) print *,'after call to getgbe, j= ',j,' k= ',k,' iret= ' & ,iret c if (iret.eq.0) then j=k call grange(kf,lb,f,dmin,dmax) print '(i4,2x,8i5,5i5,i8,2g12.4)', & k,(kpds(i),i=5,11),kpds(14),kens,kf,dmin,dmax if (convflag .eq. 'y') call conv1d2d (f,kf,data2d,nlon,nlat) if (prtpdsflag .eq. 'y' .and. iret .eq. 0) & call prtpds (kpds,kgds,k) if (newgrbflag .eq. 'y') & call gribout (lugout,kf,kpds,kgds,kens,lb,f,iret) else print *,' ' print *,'!!! ERROR in getgbe: nf=',nf,' nt= ',nt & ,' iret=',iret endif c enddo enddo c stop end c c--------------------------------------------------c c c c--------------------------------------------------c subroutine grange(n,ld,d,dmin,dmax) logical ld dimension ld(n),d(n) c dmin=1.e40 dmax=-1.e40 c do i=1,n if(ld(i)) then dmin=min(dmin,d(i)) dmax=max(dmax,d(i)) endif enddo c return end c c--------------------------------------------------c c c c--------------------------------------------------c subroutine conv1d2d (d,np,data2d,nlon,nlat) c c Convert the 1-d array (d) to a 2-d array (data2d) c real d(np),data2d(nlon,nlat) c do ilat=1,nlat do ilon=1,nlon data2d(ilon,ilat) = d(ilon+(ilat-1)*nlon) enddo enddo c return end c c---------------------------------------------------------c c c c---------------------------------------------------------c subroutine gribout (lugout,kf,kpds,kgds,kens,lb,d,iret) c c This subroutine will write out a new grib record. It does c this by calling putgbe, which is an /nwprod/w3lib routine. c integer kpds(25),kgds(22),kens(5) logical lb(kf) real d(kf) c print *,' ' print *,'at beginning of gribout, kf = ',kf c call putgbe (lugout,kf,kpds,kgds,kens,lb,d,iret) c if (iret.eq.0) then print *,' ' print *,'IRET = 0 after call to putgb, num pts written = ',kf print *,' ' else print *,' ' print *,'!!! IRET NE 0 AFTER CALL TO PUTGB, IRET= ',iret print *,' ' endif c stop end c c---------------------------------------------------------c c c c---------------------------------------------------------c subroutine prtpds (kpds,kgds,krec) c c Print out the kpds and kgds for the current record c integer kpds(25),kgds(22) c print *,' ' print *,' ==== kpds and kgds for record number ',krec,' ====' c write(*,980) kpds(1),kpds(2) write(*,981) kpds(3),kpds(4) write(*,982) kpds(5),kpds(6) write(*,983) kpds(7),kpds(8) write(*,984) kpds(9),kpds(10) write(*,985) kpds(11),kpds(12) write(*,986) kpds(13),kpds(14) write(*,987) kpds(15),kpds(16) write(*,988) kpds(17),kpds(18) write(*,989) kpds(19),kpds(20) write(*,990) kpds(21),kpds(22) write(*,991) kpds(23),kpds(24) write(*,992) kpds(25) write(*,880) kgds(1),kgds(2) write(*,881) kgds(3),kgds(4) write(*,882) kgds(5),kgds(6) write(*,883) kgds(7),kgds(8) write(*,884) kgds(9),kgds(10) write(*,885) kgds(11),kgds(12) write(*,886) kgds(13),kgds(14) write(*,887) kgds(15),kgds(16) write(*,888) kgds(17),kgds(18) write(*,889) kgds(19),kgds(20) c 980 format(' kpds(1) = ',i7,' kpds(2) = ',i7) 981 format(' kpds(3) = ',i7,' kpds(4) = ',i7) 982 format(' kpds(5) = ',i7,' kpds(6) = ',i7) 983 format(' kpds(7) = ',i7,' kpds(8) = ',i7) 984 format(' kpds(9) = ',i7,' kpds(10) = ',i7) 985 format(' kpds(11) = ',i7,' kpds(12) = ',i7) 986 format(' kpds(13) = ',i7,' kpds(14) = ',i7) 987 format(' kpds(15) = ',i7,' kpds(16) = ',i7) 988 format(' kpds(17) = ',i7,' kpds(18) = ',i7) 989 format(' kpds(19) = ',i7,' kpds(20) = ',i7) 990 format(' kpds(21) = ',i7,' kpds(22) = ',i7) 991 format(' kpds(23) = ',i7,' kpds(24) = ',i7) 992 format(' kpds(25) = ',i7) 880 format(' kgds(1) = ',i7,' kgds(2) = ',i7) 881 format(' kgds(3) = ',i7,' kgds(4) = ',i7) 882 format(' kgds(5) = ',i7,' kgds(6) = ',i7) 883 format(' kgds(7) = ',i7,' kgds(8) = ',i7) 884 format(' kgds(9) = ',i7,' kgds(10) = ',i7) 885 format(' kgds(11) = ',i7,' kgds(12) = ',i7) 886 format(' kgds(13) = ',i7,' kgds(14) = ',i7) 887 format(' kgds(15) = ',i7,' kgds(16) = ',i7) 888 format(' kgds(17) = ',i7,' kgds(18) = ',i7) 889 format(' kgds(19) = ',i7,' kgds(20) = ',i7) c return end --------------------------------------------------------------------