!> @file !! @brief Determine whether file is grib or not. !! @author gayno org: w/np2 @date 2007-nov-28 !> Determine whether file is grib or not. !! !! program history log: !! - 2007-nov-28 gayno - initial version !! - 2011-apr-26 gayno - replace my simple-minded logic !! with call to w3lib routin skgb. !! - 2014-feb-07 gayno - determine whether file is !! grib1 or grib2. !! !! @param[in] file_name - file to be checked !! @param[out] isgrib - '1' or '2' if grib1/2 file '0' if not grib !! !! input files: !! - file to be checked, fort.11 !! !! condition codes: all fatal !! - bad file open, fort.11 !! subroutine grib_check(file_name, isgrib) implicit none character*(*), intent(in) :: file_name integer, parameter :: iunit=11 integer :: istat, iseek, mseek, lskip, lgrib, version integer, intent(out) :: isgrib print*,"- CHECK FILE TYPE OF: ", trim(file_name) call baopenr (iunit, file_name, istat) if (istat /= 0) then print*,'- FATAL ERROR: BAD FILE OPEN. ISTAT IS ',istat call w3tage('SNOW2MDL') call errexit(40) end if iseek = 0 mseek = 64 call skgb2(iunit, iseek, mseek, lskip, lgrib, version) call baclose(iunit, istat) if (lgrib > 0) then isgrib = version if (isgrib == 1) print*,"- FILE IS GRIB1" if (isgrib == 2) print*,"- FILE IS GRIB2" else isgrib = 0 print*,"- FILE IS BINARY" endif return end subroutine grib_check !> Determine whether file is grib or not. !! !! Based on w3nco library routine skgb. !! !! @param[in] lugb file unit number !! @param[in] iseek number of bits to skip before search. !! @param[in] mseek max number of bytes to search. !! @param[out] lskip number of bytes to skip before message !! @param[out] lgrib number of bytes in message. '0' if not grib. !! @param[out] i1 '1' or '2' if grib1/2 file. '0' if not grib. !! !! input file: !! - file to be checked, unit=lugb !! !! @author George Gayno org: w/np2 @date 2014-Feb-07 SUBROUTINE SKGB2(LUGB,ISEEK,MSEEK,LSKIP,LGRIB,I1) implicit none INTEGER, INTENT( IN) :: LUGB, ISEEK, MSEEK INTEGER, INTENT(OUT) :: LSKIP, LGRIB, I1 INTEGER, PARAMETER :: LSEEK=128 INTEGER :: K, KZ, KS, KG, KN, KM, I4, K4 CHARACTER Z(LSEEK) CHARACTER Z4(4) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - I1=0 LGRIB=0 KS=ISEEK KN=MIN(LSEEK,MSEEK) KZ=LSEEK ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! LOOP UNTIL GRIB MESSAGE IS FOUND DO WHILE(LGRIB.EQ.0.AND.KN.GE.8.AND.KZ.EQ.LSEEK) ! READ PARTIAL SECTION CALL BAREAD(LUGB,KS,KN,KZ,Z) KM=KZ-8+1 K=0 ! LOOK FOR 'GRIB...1' IN PARTIAL SECTION DO WHILE(LGRIB.EQ.0.AND.K.LT.KM) CALL GBYTEC(Z,I4,(K+0)*8,4*8) CALL GBYTEC(Z,I1,(K+7)*8,1*8) IF(I4.EQ.1196575042.AND.(I1.EQ.1.OR.I1.EQ.2)) THEN ! LOOK FOR '7777' AT END OF GRIB MESSAGE IF (I1.EQ.1) CALL GBYTEC(Z,KG,(K+4)*8,3*8) IF (I1.EQ.2) CALL GBYTEC(Z,KG,(K+12)*8,4*8) CALL BAREAD(LUGB,KS+K+KG-4,4,K4,Z4) IF(K4.EQ.4) THEN CALL GBYTEC(Z4,I4,0,4*8) IF(I4.EQ.926365495) THEN ! GRIB MESSAGE FOUND LSKIP=KS+K LGRIB=KG ENDIF ENDIF ENDIF K=K+1 ENDDO KS=KS+KM KN=MIN(LSEEK,ISEEK+MSEEK-KS) ENDDO ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END subroutine skgb2 !> Convert from the grib2 grid description template array !! used by the ncep grib2 library, to the grib1 grid !! description section array used by ncep ipolates library. !! !! @param[in] igdtnum grib2 grid desc template number !! @param[in] igdstmpl grib2 grid desc template array !! @param[in] igdtlen grib2 grid desc template array size !! @param[out] kgds grib1 grid description section array used by ncep ipolates library. !! @param[out] ni i grid dimensions !! @param[out] nj j grid dimensions !! @param[out] res grid resolution in km !! !! condition codes: !! 50 - unrecognized model grid type; fatal !! !! @author George Gayno org: w/np2 @date 2014-Sep-26 subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) implicit none integer, intent(in ) :: igdtnum, igdtlen, igdstmpl(igdtlen) integer, intent( out) :: kgds(200), ni, nj integer :: iscale real, intent( out) :: res kgds=0 if (igdtnum.eq.0) then ! lat/lon grid iscale=igdstmpl(10)*igdstmpl(11) if (iscale == 0) iscale = 1e6 kgds(1)=0 ! oct 6 kgds(2)=igdstmpl(8) ! octs 7-8, Ni ni = kgds(2) kgds(3)=igdstmpl(9) ! octs 9-10, Nj nj = kgds(3) kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point kgds(6)=0 ! oct 17, resolution and component flags if (igdstmpl(1)==2 ) kgds(6)=64 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, di kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(19),7)) kgds(11) = 128 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32 kgds(12)=0 ! octs 29-32, reserved kgds(19)=0 ! oct 4, # vert coordinate parameters kgds(20)=255 ! oct 5, used for thinned grids, set to 255 res = float(kgds(9)) / 1000.0 * 111.0 elseif (igdtnum.eq.40) then ! Gaussian Lat/Lon grid iscale=igdstmpl(10)*igdstmpl(11) if (iscale==0) iscale=1e6 kgds(1)=4 ! oct 6 kgds(2)=igdstmpl(8) ! octs 7-8, Ni ni = kgds(2) kgds(3)=igdstmpl(9) ! octs 9-10, Nj nj = kgds(3) kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point kgds(6)=0 ! oct 17, resolution and component flags if (igdstmpl(1)==2 ) kgds(6)=64 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, Di kgds(10)=igdstmpl(18) ! octs 26-27, Number of parallels kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(19),7)) kgds(11) = 128 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32 kgds(12)=0 ! octs 29-32, reserved kgds(19)=0 ! oct 4, # vert coordinate parameters kgds(20)=255 ! oct 5, used for thinned grids, set to 255 res = float(kgds(9)) / 1000.0 * 111.0 elseif (igdtnum.eq.20) then ! Polar Stereographic Grid iscale=1e6 kgds(1)=5 ! oct 6, data representation type, polar kgds(2)=igdstmpl(8) ! octs 7-8, nx ni = kgds(2) kgds(3)=igdstmpl(9) ! octs 8-10, ny nj = kgds(3) kgds(4)=nint(float(igdstmpl(10))/float(iscale)*1000.) ! octs 11-13, lat of 1st grid point kgds(5)=nint(float(igdstmpl(11))/float(iscale)*1000.) ! octs 14-16, lon of 1st grid point kgds(6)=0 ! oct 17, resolution and component flags if (igdstmpl(1) >= 2 .or. igdstmpl(1) <= 5) kgds(6)=64 if (igdstmpl(1) == 7) kgds(6)=64 if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128 if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8 kgds(7)=nint(float(igdstmpl(14))/float(iscale)*1000.) ! octs 18-20, lon of orientation kgds(8)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 21-23, dx kgds(9)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 24-26, dy kgds(10)=0 ! oct 27, projection center flag if (btest(igdstmpl(17),1)) kgds(10) = 128 kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(18),7)) kgds(11) = 128 if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) + 64 if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) + 32 kgds(19)=0 ! oct 4, # vert coordinate parameters kgds(20)=255 ! oct 5, used for thinned grids, set to 255 res = 0.5 * float(kgds(8)+kgds(9)) / 1000. elseif (igdtnum.eq.1) then ! Rotated Lat/Lon grid if (btest(igdstmpl(19),2)) then ! e-stagger, bit 6 of scan mode is '1' iscale=igdstmpl(10)*igdstmpl(11) if (iscale == 0) iscale = 1e6 kgds(1)=203 ! oct 6, "E" grid kgds(2)=igdstmpl(8) ! octs 7-8, Ni ni = kgds(2) kgds(3)=igdstmpl(9) ! octs 9-10, Nj nj = kgds(3) kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point kgds(6)=0 ! oct 17, resolution and component flags if (igdstmpl(1)==2 ) kgds(6)=64 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8 kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000 ! octs 18-20, Lat of cent of rotation kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23, Lon of cent of rotation kgds(9)=nint(float(igdstmpl(17))/float(iscale)*500.) ! octs 24-25, Di ! Note!! grib 2 convention twice grib 1 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, Dj kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(19),7)) kgds(11) = 128 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32 kgds(12)=0 ! octs 29-32, reserved kgds(19)=0 ! oct 4, # vert coordinate parameters kgds(20)=255 ! oct 5, used for thinned grids, set to 255 res = sqrt( (float(kgds(9)) / 1000.0)**2 + & (float(kgds(10)) / 1000.0)**2 ) res = res * 111.0 else ! b-stagger iscale=igdstmpl(10)*igdstmpl(11) if (iscale == 0) iscale = 1e6 kgds(1)=205 ! oct 6, rotated lat/lon for Non-E Stagger grid kgds(2)=igdstmpl(8) ! octs 7-8, Ni ni = kgds(2) kgds(3)=igdstmpl(9) ! octs 9-10, Nj nj = kgds(3) kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point kgds(6)=0 ! oct 17, resolution and component flags if (igdstmpl(1)==2 ) kgds(6)=64 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8 kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000 ! octs 18-20, Lat of cent of rotation kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23, Lon of cent of rotation kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, Di kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, Dj kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(19),7)) kgds(11) = 128 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32 kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of last grid point kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of last grid point kgds(19)=0 ! oct 4, # vert coordinate parameters kgds(20)=255 ! oct 5, used for thinned grids, set to 255 res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) & * 0.5 * 111.0 endif else print*,'- FATAL ERROR CONVERTING TO GRIB2 GDT.' print*,'- UNRECOGNIZED GRID TYPE.' call w3tage('SNOW2MDL') call errexit(50) endif end subroutine gdt_to_gds !> Determine length of grib2 gds template array, which is a function of !! the map projection. !! !! @note call this routine before init_grib2. !! !! @param[in] kgds grib1 gds array !! @param[in] igdstmplen length of gds template array. !! !! condition codes: !! 47 - unrecognized grid type; fatal !! !! @author George Gayno org: w/np2 @date 2014-Sep-28 subroutine grib2_check (kgds, igdstmplen) implicit none integer, intent(in) :: kgds(200) integer, intent( out) :: igdstmplen select case (kgds(1)) case(4) ! gaussian igdstmplen = 19 case(203, 205) ! rotated lat/lon "B" or "E" stagger igdstmplen = 22 case default print*,'- FATAL ERROR IN ROUTINE GRIB2_CHECK.' print*,'- UNRECOGNIZED GRID TYPE.' call w3tage('SNOW2MDL') call errexit(47) end select end subroutine grib2_check !> Initialize grib2 arrays required by the ncep g2 library according to !! grib1 gds information. The grib1 gds is held in the kgds array, which !! is used by the ncep ipolates and w3nco (grib 1) libraries. !! !! Call routine grib2_check first to determine igdstmplen. !! !! @param[in] century current date/time info !! @param[in] year current date/time info !! @param[in] month current date/time info !! @param[in] day current date/time info !! @param[in] hour current date/time info !! @param[in] kgds grib1 gds information !! @param[in] igdstmplen length of grib2 gdt template. !! @param[in] lat11 lat of first grid point !! @param[in] lon11 lon of first grid point !! @param[in] latlast lat of last grid point !! @param[in] lonlast lon of last grid point !! @param[out] igds grib2 section 3 information. !! @param[out] listsec0 grib2 section 0 information. !! @param[out] listsec1 grib2 section 1 information. !! @param[out] ipdsnum grib2 pds template number !! @param[out] ipdstmpl grib2 pds template array !! @param[out] igdstmpl grib2 gds template array !! @param[out] idefnum information for non-reg grid, grid points in each row. !! @param[out] ideflist information for non-reg grid, grid points in each row. !! @param[out] ngrdpts number of model grid points. !! @author George Gayno org: w/np2 @date 2014-Sep-28 subroutine init_grib2(century, year, month, day, hour, kgds, & lat11, latlast, lon11, lonlast, & listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl, & igdstmplen, idefnum, ideflist, ngrdpts) implicit none integer, intent(in ) :: century, year, month, day, hour integer, intent(in ) :: kgds(200), igdstmplen integer, intent( out) :: igds(5) integer, intent( out) :: listsec0(2) integer, intent( out) :: listsec1(13) integer, intent( out) :: ipdstmpl(15), ipdsnum integer, intent( out) :: igdstmpl(igdstmplen) integer, intent( out) :: idefnum, ideflist integer, intent( out) :: ngrdpts real, intent(in ) :: lat11, latlast, lon11, lonlast real :: scale ! Section 0 listsec0(1)=0 ! discipline, meteorological fields listsec0(2)=2 ! grib version 2 ! Section 1 listsec1(1)=7 ! id of center (ncep) listsec1(2)=4 ! subcenter (emc) listsec1(3)=8 ! master table version number. wgrib2 does not recognize later tables listsec1(4)= 0 ! local table not used listsec1(5)= 0 ! signif of ref time - analysis if (year == 100) then listsec1(6)=century*100 + year else listsec1(6)=(century-1)*100 + year endif listsec1(7)=month listsec1(8)=day listsec1(9)=hour listsec1(10:11)=0 ! minutes/secs listsec1(12)=0 ! production status of data - ops products listsec1(13)=0 ! type of processed products - analysis ! Section 2 - not used ! Section 3 - grid description section if (kgds(1) == 4) then ! gaussian igdstmpl(1)=5 ! oct 15; shape of the earth, wgs84 igdstmpl(2)=255 ! oct 16; scale factor of radius of spherical earth, not used. igdstmpl(3)=-1 ! octs 17-20; scale value of radius of spherical earth, not used. igdstmpl(4)=255 ! oct 21; scale factor of major axis of elliptical earth, not used. igdstmpl(5)=-1 ! octs 22-25; scaled value of major axis of elliptical earth, not used. igdstmpl(6)=255 ! oct 26; scale factor of minor axis of elliptical earth, not used. igdstmpl(7)=-1 ! octs 27-30; scaled value of minor axis of elliptical earth, not used. igdstmpl(8)=kgds(2) ! octs 31-34; # "i" points igdstmpl(9)=kgds(3) ! octs 35-38; # "j" points igdstmpl(10)=1 ! octs 39-42; basic angle igdstmpl(11)=10**6 ! octs 43-46; subdivisions of basic angle scale=float(igdstmpl(10)*igdstmpl(11)) igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point if (lon11 < 0) then igdstmpl(13)=nint((lon11+360.)*scale) ! octs 51-54; lon of first grid point else igdstmpl(13)=nint(lon11*scale) endif igdstmpl(14) = 0 ! oct 55; resolution and component flags if (btest(kgds(6),7)) igdstmpl(14) = 48 if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8 igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point if (lonlast < 0) then igdstmpl(16)=nint((lonlast+360.)*scale) ! octs 60-63; lon of last grid point else igdstmpl(16)=nint(lonlast*scale) endif igdstmpl(17)= nint(360.0/float(kgds(2)-1)*scale) ! octs 64-67; di of grid igdstmpl(18)= kgds(3)/2 ! octs 68-71; # grid pts between pole and equator igdstmpl(19)=0 ! oct 72; scanning mode flag if(btest(kgds(11),7)) igdstmpl(19)=128 if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64 if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32 igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1 igds(2) = kgds(2)*kgds(3) ! num grid points igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid) igds(4) = 0 ! regular grid, no appended list igds(5) = 40 ! gaussian ngrdpts = igds(2) ! These variables used for non-regular grids. We are using regular grids ! (igds(3) equals 0). idefnum=1 ideflist=0 elseif (kgds(1) == 203 .or. kgds(1) == 205) then igdstmpl(1)=5 ! oct 15; shape of the earth, wgs84 igdstmpl(2)=255 ! oct 16; scale factor of radius of spherical earth, not used. igdstmpl(3)=-1 ! octs 17-20; scale value of radius of spherical earth, not used. igdstmpl(4)=255 ! oct 21; scale factor of major axis of elliptical earth, not used. igdstmpl(5)=-1 ! octs 22-25; scaled value of major axis of elliptical earth, not used. igdstmpl(6)=255 ! oct 26; scale factor of minor axis of elliptical earth, not used. igdstmpl(7)=-1 ! octs 27-30; scaled value of minor axis of elliptical earth, not used. igdstmpl(8)=kgds(2) ! octs 31-34; # "i" points igdstmpl(9)=kgds(3) ! octs 35-38; # "j" points igdstmpl(10)=1 ! octs 39-42; basic angle igdstmpl(11)=10**6 ! octs 43-46; subdivisions of basic angle scale=float(igdstmpl(10)*igdstmpl(11)) igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point if (lon11 < 0) then igdstmpl(13)=nint((lon11+360.)*scale) ! octs 51-54; lon of first grid point else igdstmpl(13)=nint(lon11*scale) endif igdstmpl(14) = 0 ! oct 55; resolution and component flags if (btest(kgds(6),7)) igdstmpl(14) = 48 if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8 igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point if (lonlast < 0) then igdstmpl(16)=nint((lonlast+360.)*scale) ! octs 60-63; lon of last grid point else igdstmpl(16)=nint(lonlast*scale) endif if (kgds(1) == 203) igdstmpl(17)= nint(float(kgds(9))*scale/500.) ! octs 64-67; di of grid. ! iplib "e" grid convention ! is 1/2 the grib convention. if (kgds(1) == 205) igdstmpl(17)= nint(float(kgds(9))*scale/1000.) ! octs 64-67; di of grid igdstmpl(18)= nint(float(kgds(10))*scale/1000.) ! octs 68-71; dj of grid igdstmpl(19)=0 ! oct 72; scanning mode flag if(btest(kgds(11),7)) igdstmpl(19)=128 if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64 if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32 if (kgds(1) == 203) igdstmpl(19)=igdstmpl(19) + 4 igdstmpl(20) = nint(float(kgds(7)-90000)*scale/1000.) ! octs 73-76; lat of south pole of projection if (kgds(8) < 0) then igdstmpl(21) = nint(float(kgds(8)+360000)*scale/1000.) ! octs 77-80; long of southern pole of projection. else igdstmpl(21) = nint(float(kgds(8))*scale/1000.) ! octs 77-80; long of southern pole of projection. endif igdstmpl(22)=0 ! octs 81-84; angle of rotation of projection igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1 igds(2) = kgds(2)*kgds(3) ! num grid points igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid) igds(4) = 0 ! regular grid, no appended list igds(5) = 1 ! rotated lat/lon ngrdpts = igds(2) ! These variables used for non-regular grids. We are using regular grids ! (igds(3) equals 0). idefnum=1 ideflist=0 end if ! Section 4 - product definition section ipdsnum = 0 ! pds template number - table 4.0 ipdstmpl(1)= 1 ! oct 10; parameter category ! note!! to use a parmeter number >= 192 you must set the local table to '1' ipdstmpl(2)= 42 ! oct 11; parameter ipdstmpl(3)= 0 ! oct 12; type of generating process ipdstmpl(4)= 255 ! oct 13; background generating process identifier ipdstmpl(5)= 84 ! oct 14; analysis generating process identifier ipdstmpl(6)= 0 ! octs 15-16; hours after ob cutoff ipdstmpl(7)= 0 ! oct 17; minutes after ob cutoff ipdstmpl(8)= 1 ! oct 18; unit of time range ipdstmpl(9)= 0 ! octs 19-22; forecast time in units defined by oct 18 ipdstmpl(10)=1 ! oct 23; type of first fixed surface ipdstmpl(11)=0 ! oct 24; scale factor of first fixed surface ipdstmpl(12)=0 ! octs 25-28; scale value of first fixed surface ipdstmpl(13)=255 ! oct 29; type of second fixed surface ipdstmpl(14)=255 ! oct 30; scale factor of second fixed surface ipdstmpl(15)=-2147483647 ! octs 31-34; scaled value of second fixed surface ! note! for these particular octets, using -1 as ! missing does not work because -1 may be an actual ! scaled value. after looking thru the g2 library ! and some trial and error, i determined that missing ! is minus 2**31-1. end subroutine init_grib2 !> Nullify the grib2 gribfield pointers. !! !! @param[in] gfld a gribfield data structure !! !! @author George Gayno org: w/np2 @date 2014-Sep-28 subroutine grib2_null(gfld) use grib_mod implicit none type(gribfield), intent(inout) :: gfld nullify(gfld%idsect) nullify(gfld%local) nullify(gfld%list_opt) nullify(gfld%igdtmpl) nullify(gfld%ipdtmpl) nullify(gfld%coord_list) nullify(gfld%idrtmpl) nullify(gfld%bmap) nullify(gfld%fld) end subroutine grib2_null !> Deallocate the grib2 gribfield pointers. !! !! @param[in] gfld a gribfield data structure !! !! @author George Gayno org: w/np2 @date 2014-Sep-28 subroutine grib2_free(gfld) use grib_mod implicit none type(gribfield), intent(inout) :: gfld if (associated(gfld%idsect)) deallocate(gfld%idsect) if (associated(gfld%local)) deallocate(gfld%local) if (associated(gfld%list_opt)) deallocate(gfld%list_opt) if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl) if (associated(gfld%ipdtmpl)) deallocate(gfld%ipdtmpl) if (associated(gfld%coord_list)) deallocate(gfld%coord_list) if (associated(gfld%idrtmpl)) deallocate(gfld%idrtmpl) if (associated(gfld%bmap)) deallocate(gfld%bmap) if (associated(gfld%fld)) deallocate(gfld%fld) end subroutine grib2_free