! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE RDNMCGRB2 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE rdnmcgrb2(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, & 4,26 iproj_grb,nx_grb,ny_grb,dx_grb,dy_grb, & latsw,lonsw,lattru1,lattru2,lontrue,uvearth, & n2dvs, n3dvs, maxvar, nzsoilin_ext, & varids, var2dindx, var3dindx, var2dlvl, var3dlvl, var3dsoil, & var_grb2d, var_grb3d, lvldbg, & iret) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! This subroutine is to read and select given variables with GRIB ! grid ID and level type from NMC GRIB files. ! ! The decoder of GRIB2 is from NCEP. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 08/01/2006 ! ! MODIFICATIONS: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! gribfile GRIB file name ! grbflen Length of GRIB file name ! gribtime Time of GRIB data ! ! OUTPUT: ! ! iproj_grb Map projection number of GRIB data ! trlon_grb True longitude of GRIB data (degrees E) ! latnot_grb(2) True latitude(s) of GRIB data (degrees N) ! swlat_grb Latitude of first grid point at southwest conner ! swlon_grb Longitude of first grid point at southwest conner ! ! dx_grb x-direction grid length ! dy_grb y-direction grid length ! ! iret Return flag ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! USE GRIB_MOD IMPLICIT NONE INCLUDE 'mp.inc' INTEGER, INTENT(IN) :: nx_ext,ny_ext,nz_ext INTEGER, INTENT(IN) :: grbflen CHARACTER(LEN=*), INTENT(IN) :: gribfile CHARACTER(LEN=*), INTENT(IN) :: gribtime INTEGER, INTENT(OUT) :: iproj_grb ! Map projection indicator ! Already converted to ARPS map definitions INTEGER, INTENT(OUT) :: nx_grb ! Number of points along x-axis INTEGER, INTENT(OUT) :: ny_grb ! Number of points along y-axis REAL, INTENT(OUT) :: dx_grb ! x-direction increment or grid length REAL, INTENT(OUT) :: dy_grb ! y-direction increment or grid length REAL, INTENT(OUT) :: latsw ! Latitude of South West corner point REAL, INTENT(OUT) :: lonsw ! Longitude of South West corner point REAL, INTENT(OUT) :: lattru1 ! Latitude (1st) at which projection is true REAL, INTENT(OUT) :: lattru2 ! Latitude (2nd) at which projection is true REAL, INTENT(OUT) :: lontrue ! Longitude at which projection is true INTEGER, INTENT(OUT) :: uvearth ! = 0, Resolved u and v components of vector quantities relative to easterly and northerly directions ! = 1, Resolved u and v components of vector quantities relative to the defined grid in the direction of increasing x and y (or i and j) coordinates, respectively INTEGER, INTENT(IN) :: n2dvs, n3dvs, maxvar INTEGER, INTENT(IN) :: nzsoilin_ext INTEGER, INTENT(IN) :: varids(4,maxvar) INTEGER, INTENT(IN) :: var2dindx(maxvar), var3dindx(maxvar) REAL, INTENT(IN) :: var2dlvl(n2dvs), var3dlvl(nz_ext), var3dsoil(nzsoilin_ext) REAL, INTENT(OUT) :: var_grb2d(nx_ext,ny_ext,n2dvs) REAL, INTENT(OUT) :: var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs) INTEGER, INTENT(IN) :: lvldbg INTEGER, INTENT(OUT) :: iret ! Return flag ! !----------------------------------------------------------------------- ! ! Temporary arrays to read GRIB file ! !----------------------------------------------------------------------- ! INTEGER :: grbunit INTEGER :: istatus ! Return status INTEGER, PARAMETER :: maxlen = 32000 CHARACTER(LEN=1), ALLOCATABLE :: cgrib(:) TYPE(gribfield) :: gfld ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,m,n,ij INTEGER :: iyr, imo, iday, ihr, fhr LOGICAL :: fexist INTEGER :: grbflen_new CHARACTER (LEN=256) :: gribfile_new LOGICAL :: verbose = .FALSE. CHARACTER(LEN=12), PARAMETER :: cmd = 'rdnmcgrib2: ' INTEGER :: iseek, icount, itot INTEGER :: lgrib, lgribin, lskip INTEGER :: currlen INTEGER :: listsec0(3), listsec1(13) INTEGER :: numfields, numlocal, maxlocal REAL :: levelin INTEGER :: iscan, jscan, ijscan INTEGER :: ibgn, iinc, iend, jbgn, jinc, jend INTEGER :: itemp ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF (lvldbg > 90) verbose = .TRUE. !----------------------------------------------------------------------- ! ! Date and time expected ! !----------------------------------------------------------------------- READ (gribtime,'(i4,3i2,1x,i3)') iyr,imo,iday,ihr,fhr ! !----------------------------------------------------------------------- ! ! Open the GRIB file ! !----------------------------------------------------------------------- ! INQUIRE(FILE=gribfile(1:grbflen), EXIST = fexist ) IF( fexist ) THEN gribfile_new = gribfile grbflen_new = grbflen GO TO 200 ENDIF INQUIRE(FILE=trim(gribfile(1:grbflen))//'.Z', EXIST = fexist ) IF( fexist ) THEN CALL unixcmd( 'cp '//trim(gribfile(1:grbflen))//'.Z tem_ext_file.Z' ) INQUIRE(FILE='tem_ext_file', EXIST = fexist ) IF( fexist ) call unixcmd( '/bin/rm tem_ext_file' ) CALL uncmprs( 'tem_ext_file.Z' ) gribfile_new = 'tem_ext_file' grbflen_new = 12 GO TO 200 END IF INQUIRE(FILE=trim(gribfile(1:grbflen))//'.gz', EXIST = fexist ) IF( fexist ) THEN CALL unixcmd( 'cp '//trim(gribfile(1:grbflen))//'.gz tem_ext_file.gz' ) INQUIRE(FILE='tem_ext_file', EXIST = fexist ) IF( fexist ) call unixcmd( '/bin/rm tem_ext_file' ) CALL uncmprs( 'tem_ext_file.gz' ) gribfile_new = 'tem_ext_file' grbflen_new = 12 GO TO 200 END IF WRITE(6,'(/1x,a,/1x,a/)') WRITE(6,'(a)') 'GRIB file '//gribfile(1:grbflen)// & ' or its compressed version not found. Program stopped in RDNMCGRB2.' STOP 200 CONTINUE WRITE(6,'(1x,8A)') cmd, 'Opening GRIB ', gribfile(1:grbflen) CALL getunit( grbunit ) CALL baopenr(grbunit,gribfile_new(1:grbflen_new),iret) IF ( iret /= 0 ) THEN WRITE(6,'(a,a,/,a,i4)') & 'ERROR in opening file ',gribfile_new(1:grbflen_new), & 'function baopenr return status: ',iret RETURN END IF !----------------------------------------------------------------------- ! ! Loop over all records ! !----------------------------------------------------------------------- itot = 0 ! total field counter icount = 0 ! total GRIB2 message counter iseek = 0 currlen = 0 ! DO n = 1,maxvar ! WRITE(0,*) varids(:,n) ! END DO ! STOP iret = 0 DO WHILE (iret == 0) CALL skgb(grbunit,iseek,maxlen,lskip,lgrib) IF (lgrib == 0) EXIT IF (lgrib > currlen) THEN ! Extend buffer size as needed IF (ALLOCATED(cgrib)) DEALLOCATE(cgrib) ALLOCATE(cgrib(lgrib), STAT = istatus) currlen = lgrib END IF CALL baread(grbunit,lskip,lgrib,lgribin,cgrib) ! read a GRIB message IF (lgrib /= lgribin) THEN WRITE(6,'(2(a,I6))') & 'ERROR: GRIB2 message size in error, required size = ',lgrib, & ', read-in size = ',lgribin CALL arpsstop('ERROR in rdnmcgrb2 when calling baread.',1) END IF icount = icount + 1 iseek = lskip+lgrib ! next GRIB message starting location CALL gb_info(cgrib,lgrib,listsec0,listsec1, & numfields,numlocal,maxlocal,iret) IF (iret /= 0) THEN WRITE(6,'(a,I3)') 'ERROR: quering GRIB2 message = ',iret CALL arpsstop('ERROR in rdnmcgrb2 when calling gb_info.',1) END IF itot = itot+numfields IF (icount == 1) THEN ! Once per file for date and map projection value CALL gf_getfld(cgrib,lgrib,1,.TRUE.,.TRUE.,gfld,iret) IF (iret /= 0) THEN WRITE(6,'(a,I3)') 'ERROR: Decoding GRIB2 message = ',iret CALL arpsstop('ERROR in rdnmcgrb2 when calling gf_getfld.',1) END IF IF (verbose) & WRITE(6,'(1x,a,I4.4,4(a,I2.2),/)') 'GRIB date-time found is: ', & gfld%idsect(6),'-',gfld%idsect(7),'-',gfld%idsect(8),'_', & gfld%idsect(9),'f',gfld%ipdtmpl(9) IF (iyr /= gfld%idsect(6) .OR. imo /= gfld%idsect(7) .OR. & iday /= gfld%idsect(8).OR. ihr /= gfld%idsect(9) .OR. & fhr /= gfld%ipdtmpl(9) ) THEN WRITE(6,'(a,/,7x,2a,/,7x,a,I4.4,4(a,I2.2))') & 'ERROR: GRIB data is in wrong date and time', & 'Expected: ',gribtime, ' Found: ',gfld%idsect(6),'-', & gfld%idsect(7),'-',gfld%idsect(8),'_',gfld%idsect(9), & 'f',gfld%ipdtmpl(9) CALL arpsstop('Wrong grib file?',1) END IF itemp = gfld%ipdtmpl(5) ! indicator of model IF (gfld%idsect(1) == 7) THEN IF (itemp == 83 .OR. itemp == 84) THEN WRITE(6,'(1x,a,/)') 'NCEP NAM model' ELSE IF (itemp == 81 .OR. itemp == 96) THEN WRITE(6,'(1x,a,/)') 'NCEP GFS model' ELSE WRITE(6,'(1x,a,I3,/)') 'Unkown model from NCEP, model = ',itemp END IF ELSE WRITE(6,'(1x,a,I2,a,I2,/)') & 'Unknown model and orig center, center = ',gfld%idsect(1), & ' model = ',itemp END IF IF (gfld%igdtnum == 0) THEN ! Lat/Lon grid iproj_grb = 4 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) dx_grb = gfld%igdtmpl(17)/1.0E6 ! in degree dy_grb = gfld%igdtmpl(18)/1.0E6 ! in degree iscan = IAND(gfld%igdtmpl(19),128)/128 ! get the first bit ! = 0, +x(+i) direction, = 1, -x(-i) direction jscan = IAND(gfld%igdtmpl(19),64) / 64 ! get the second bit ! = 0, -y(-j) direction, = 1, +y(+j) direction ijscan = IAND(gfld%igdtmpl(19),32) / 32 ! get the third bit ! = 0, x direction first, = 1, y direction first uvearth = IAND(gfld%igdtmpl(14),8) / 8 ! get the fifth bit, big-endian IF (jscan == 1) THEN ! +y(+j) direction latsw = gfld%igdtmpl(12)/1.0E6 lonsw = gfld%igdtmpl(13)/1.0E6 ELSE ! -y(-j) direction latsw = gfld%igdtmpl(15)/1.0E6 lonsw = gfld%igdtmpl(13)/1.0E6 END IF ELSE IF (gfld%igdtnum == 30) THEN ! Lambert Conformal Grid iproj_grb = 2 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) lontrue = gfld%igdtmpl(14)/1.0E6 lattru1 = gfld%igdtmpl(19)/1.0E6 lattru2 = gfld%igdtmpl(20)/1.0E6 dx_grb = gfld%igdtmpl(15)/1.0E3 dy_grb = gfld%igdtmpl(16)/1.0E3 latsw = gfld%igdtmpl(10)/1.0E6 lonsw = gfld%igdtmpl(11)/1.0E6 iscan = IAND(gfld%igdtmpl(18),128)/128 ! get the first bit ! = 0, +x(+i) direction, = 1, -x(-i) direction jscan = IAND(gfld%igdtmpl(18),64) / 64 ! get the second bit ! = 0, -y(-j) direction, = 1, +y(+j) direction ijscan = IAND(gfld%igdtmpl(18),32) / 32 ! get the third bit ! = 0, x direction first, = 1, y direction first ELSE IF (gfld%igdtnum == 20) THEN ! Polar-Stereographic Grid iproj_grb = 1 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) lattru1 = 60. lattru2 = 91. latsw = gfld%igdtmpl(10) lonsw = gfld%igdtmpl(11) ELSE WRITE(6,'(1x,a,I3)') 'Unkown projection: ',gfld%igdtnum CALL arpsstop('Unknown map project in GRIB2 file.',1) END IF iinc = (-1)**iscan ibgn = (nx_ext)**iscan iend = ibgn + (nx_ext-1)*iinc jinc = (-1)*(-1)**jscan jend = (ny_ext)**jscan jbgn = jend - (ny_ext-1)*jinc IF (nx_ext /= nx_grb .OR. ny_ext /= ny_grb) THEN WRITE(6,'(1x,a,/,8x,a,I3,a,I3,/,8x,a,I3,a,I3)') & 'ERROR: Dimension size inconsistent', & 'Expected: nx_ext = ',nx_ext,', ny_ext = ',ny_ext, & 'Found : nx_grb = ',nx_grb,', ny_grb = ',ny_grb CALL arpsstop('Inconsistent dimension size in GRIB2 file.',1) END IF IF (verbose) THEN WRITE(6,'(1x,2a)') 'M_No StartBytes F_No. ', & 'Field_No. (Dis, Cat, Par) LevelType Levelval GRB_Array' WRITE(6,'(1x,2a)') '===== ========== ===== ', & '========== =============== ========== ========== ==========' END IF CALL gf_free( gfld ) END IF ! First GRIB2 message IF (verbose) WRITE(6,FMT='(1x,I4,2x,I10,1x,I3,3x)',ADVANCE='NO') & icount,lskip+1,numfields DO n = 1,numfields ! unpack GRIB2 fields in each message CALL gf_getfld(cgrib,lgrib,n,.TRUE.,.TRUE.,gfld,iret) IF (iret /= 0) THEN WRITE(6,'(a,I3)') 'ERROR: Decoding GRIB2 message = ',iret CALL arpsstop('ERROR in rdnmcgrb2 when calling gf_getfld.',1) END IF IF (verbose) THEN IF (n == 1) THEN WRITE(6,FMT='(I5,6x,3(a,I3),a,I10,4x)',ADVANCE='NO') n, & '(',gfld%discipline,',',gfld%ipdtmpl(1),',',gfld%ipdtmpl(2),')',& gfld%ipdtmpl(10) ELSE WRITE(6,FMT='(24x,I5,6x,3(a,I3),a,I10,4x)',ADVANCE='NO') n, & '(',gfld%discipline,',',gfld%ipdtmpl(1),',',gfld%ipdtmpl(2),')',& gfld%ipdtmpl(10) END IF END IF DO m = 1,maxvar ! match the required variables ! IF (gfld%ipdtmpl(4) == varids(1,m) .AND. & ! Discipline IF (gfld%discipline == varids(1,m) .AND. & ! Discipline gfld%ipdtmpl(1) == varids(2,m) .AND. & ! Category gfld%ipdtmpl(2) == varids(3,m) .AND. & ! Parameter gfld%ipdtmpl(10)== varids(4,m) ) THEN ! Elevation levelin = gfld%ipdtmpl(12) !----------------------------------------------------------------------- ! ! Check 3D levels ! !----------------------------------------------------------------------- IF (gfld%ipdtmpl(10) == 100 .OR. & ! Pressure level gfld%ipdtmpl(10) == 105 & ! Hybrid level ! gfld%ipdtmpl(10) == 104 & ! Sigma level ) THEN ! 3D levels DO k = 1,nz_ext IF (levelin == var3dlvl(k)) EXIT END DO IF (k > nz_ext) THEN WRITE(6,'(/,1x,a,F15.2,a,/,1x,a,/)') & 'WARNING: varaible level (',levelin,') is not found in var3dlvl.', & ' Please check the parameters array passing in to rdnmcgrb2.' EXIT ! CALL arpsstop('Unknown variable level.',1) ELSE IF (verbose) WRITE(6,FMT='(F10.2,5x,a,I3)',ADVANCE='NO') & levelin, '3D-',k END IF ij = 0 DO j = jbgn, jend, jinc DO i = ibgn, iend, iinc ij = ij+1 ! ij = j + (i-1)*ny_ext ! ij = i + (j-1)*nx_ext var_grb3d(i,j,k,var3dindx(m)) = gfld%fld(ij) ! store the unpacked 2D slab END DO END DO !----------------------------------------------------------------------- ! ! Check Soil levels ! !----------------------------------------------------------------------- ELSE IF (gfld%ipdtmpl(10) == 106 ) THEN ! Depth below land surface DO k = 1,nzsoilin_ext IF (levelin == var3dsoil(k)) EXIT END DO IF (k > nzsoilin_ext) THEN WRITE(6,'(/,1x,a,F15.2,a,/,1x,a,/)') & 'WARNING: varaible level (',levelin,') is not found in var3dsoil.', & ' Please check the parameters array passing in to rdnmcgrb2.' EXIT ! CALL arpsstop('Unknown soil variable level.',1) ELSE IF (verbose) WRITE(6,FMT='(F10.2,5x,a,I3.3)',ADVANCE='NO') & levelin, '3DSOIL-',k END IF ij = 0 ! assume ijscan = 0 DO j = jbgn, jend, jinc DO i = ibgn, iend, iinc ij = ij+1 var_grb3d(i,j,k,var3dindx(m)) = gfld%fld(ij) ! store the unpacked 2D slab END DO END DO !----------------------------------------------------------------------- ! ! Check 2D level, Should be at specified 2D slab ! !----------------------------------------------------------------------- ELSE IF (gfld%ipdtmpl(10) == 103 .OR. & ! Specified height above ground gfld%ipdtmpl(10) == 104 ) THEN ! Sigma level IF (levelin == var2dlvl(var2dindx(m)) ) THEN ij = 0 DO j = jbgn, jend, jinc DO i = ibgn, iend, iinc ij = ij+1 var_grb2d(i,j,var2dindx(m)) = gfld%fld(ij) END DO END DO IF (verbose) WRITE(6,FMT='(F10.2,5x,a)',ADVANCE='NO') & levelin, '2D' ELSE IF (verbose) WRITE(6,FMT='(F10.2,5x,a)',ADVANCE='NO') & levelin, '-- Skipped' END IF ! 2D variable at the right layer !----------------------------------------------------------------------- ! ! Check 2D level, Do not need to check the level values ! !----------------------------------------------------------------------- ELSE IF (gfld%ipdtmpl(10) == 1 .OR. & ! Ground or water surface gfld%ipdtmpl(10) == 101 ) THEN ! Mean Sea Level ij = 0 DO j = jbgn, jend, jinc DO i = ibgn, iend, iinc ij = ij+1 var_grb2d(i,j,var2dindx(m)) = gfld%fld(ij) END DO END DO IF (verbose) WRITE(6,FMT='(F10.2,5x,a,I2,a)',ADVANCE='NO') & levelin, '2D-ground(',var2dindx(m),')' ELSE IF (verbose) WRITE(6,FMT='(F10.2,5x,a)',ADVANCE='NO') & levelin, '-- Skipped' END IF ! Store if matched levels also EXIT ! already matched, so exit from the required variable loop END IF ! matched var. ids END DO ! matching loop all required variables IF (verbose ) THEN IF ( m > maxvar) THEN WRITE(6,'(10x,a)') ' --- skipped.' ELSE WRITE(6,*) END IF END IF CALL gf_free( gfld ) END DO ! numfields END DO ! number of GRIB2 messages !----------------------------------------------------------------------- ! ! Close file and release the gribfield before returning ! !----------------------------------------------------------------------- CALL baclose( grbunit,iret ) CALL retunit( grbunit ) DEALLOCATE(cgrib) RETURN END SUBROUTINE rdnmcgrb2 ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE RDGRB2DIMS ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE rdgrb2dims(gribfile,grbflen,nx_grb,ny_grb,iret) 1,11 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! This subroutine extracts dimension sizes from GRIB2 file ! ! The decoder of GRIB2 is from NCEP. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 08/20/2007 ! ! MODIFICATIONS: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! gribfile GRIB file name ! grbflen Length of GRIB file name ! ! OUTPUT: ! ! nx_grb x-direction grid length ! ny_grb y-direction grid length ! ! iret Return flag ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! USE GRIB_MOD IMPLICIT NONE INCLUDE 'mp.inc' INTEGER, INTENT(IN) :: grbflen CHARACTER(LEN=*), INTENT(IN) :: gribfile INTEGER, INTENT(OUT) :: nx_grb ! Number of points along x-axis INTEGER, INTENT(OUT) :: ny_grb ! Number of points along y-axis INTEGER, INTENT(OUT) :: iret ! Return flag ! !----------------------------------------------------------------------- ! ! Temporary arrays to read GRIB file ! !----------------------------------------------------------------------- ! INTEGER :: grbunit INTEGER :: istatus ! Return status INTEGER, PARAMETER :: maxlen = 32000 CHARACTER(LEN=1), ALLOCATABLE :: cgrib(:) TYPE(gribfield) :: gfld ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! ! LOGICAL :: verbose = .FALSE. CHARACTER(LEN=12), PARAMETER :: cmd = 'rdgrib2dims: ' INTEGER :: iseek INTEGER :: lgrib, lgribin, lskip INTEGER :: iscan, jscan, ijscan INTEGER :: listsec0(3), listsec1(13) INTEGER :: numfields, numlocal, maxlocal INTEGER :: itemp ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! ! Open the GRIB file ! !----------------------------------------------------------------------- ! WRITE(6,'(1x,8A)') cmd, 'Opening GRIB ', gribfile(1:grbflen) CALL getunit( grbunit ) CALL baopenr(grbunit,gribfile(1:grbflen),iret) IF ( iret /= 0 ) THEN WRITE(6,'(a,a,/,a,i4)') & 'ERROR in opening file ',gribfile(1:grbflen), & 'function baopenr return status: ',iret RETURN END IF !----------------------------------------------------------------------- ! ! Loop over all records ! !----------------------------------------------------------------------- iseek = 0 iret = 0 CALL skgb(grbunit,iseek,maxlen,lskip,lgrib) IF (lgrib == 0) THEN iret = -3 RETURN END IF ALLOCATE(cgrib(lgrib), STAT = istatus) CALL baread(grbunit,lskip,lgrib,lgribin,cgrib) ! read a GRIB message IF (lgrib /= lgribin) THEN WRITE(6,'(2(a,I6))') & 'ERROR: GRIB2 message size in error, required size = ',lgrib, & ', read-in size = ',lgribin CALL arpsstop('ERROR in rdgrb2dims when calling baread.',1) END IF CALL gb_info(cgrib,lgrib,listsec0,listsec1, & numfields,numlocal,maxlocal,iret) IF (iret /= 0) THEN WRITE(6,'(a,I3)') 'ERROR: quering GRIB2 message = ',iret CALL arpsstop('ERROR in rdgrb2dims when calling gb_info.',1) END IF CALL gf_getfld(cgrib,lgrib,1,.TRUE.,.TRUE.,gfld,iret) IF (iret /= 0) THEN WRITE(6,'(a,I3)') 'ERROR: Decoding GRIB2 message = ',iret CALL arpsstop('ERROR in rdgrb2dims when calling gf_getfld.',1) END IF itemp = gfld%ipdtmpl(5) ! indicator of model IF (gfld%idsect(1) == 7) THEN IF (itemp == 83 .OR. itemp == 84) THEN WRITE(6,'(1x,a,/)') 'NCEP NAM model' ELSE IF (itemp == 81 .OR. itemp == 96) THEN WRITE(6,'(1x,a,/)') 'NCEP GFS model' ELSE WRITE(6,'(1x,a,I3,/)') 'Unkown model from NCEP, model = ',itemp END IF ELSE WRITE(6,'(1x,a,I2,a,I2,/)') & 'Unknown model and orig center, center = ',gfld%idsect(1), & ' model = ',itemp END IF IF (gfld%igdtnum == 0) THEN ! Lat/Lon grid ! iproj_grb = 4 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) ! dx_grb = gfld%igdtmpl(17)/1.0E6 ! in degree ! dy_grb = gfld%igdtmpl(18)/1.0E6 ! in degree ! iscan = IAND(gfld%igdtmpl(19),128)/128 ! get the first bit ! = 0, +x(+i) direction, = 1, -x(-i) direction ! jscan = IAND(gfld%igdtmpl(19),64) / 64 ! get the second bit ! = 0, -y(-j) direction, = 1, +y(+j) direction ! ijscan = IAND(gfld%igdtmpl(19),32) / 32 ! get the third bit ! = 0, x direction first, = 1, y direction first ! uvearth = IAND(gfld%igdtmpl(14),8) / 8 ! get the fifth bit, big-endian ! IF (jscan == 1) THEN ! +y(+j) direction ! latsw = gfld%igdtmpl(12)/1.0E6 ! lonsw = gfld%igdtmpl(13)/1.0E6 ! ELSE ! -y(-j) direction ! latsw = gfld%igdtmpl(15)/1.0E6 ! lonsw = gfld%igdtmpl(13)/1.0E6 ! END IF ELSE IF (gfld%igdtnum == 30) THEN ! Lambert Conformal Grid ! iproj_grb = 2 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) ! lontrue = gfld%igdtmpl(14) ! lattru1 = gfld%igdtmpl(19) ! lattru2 = gfld%igdtmpl(20) ! dx_grb = gfld%igdtmpl(15) ! dy_grb = gfld%igdtmpl(16) ! latsw = gfld%igdtmpl(10) ! lonsw = gfld%igdtmpl(11) ELSE IF (gfld%igdtnum == 20) THEN ! Polar-Stereographic Grid ! iproj_grb = 1 nx_grb = gfld%igdtmpl(8) ny_grb = gfld%igdtmpl(9) ! lattru1 = 60. ! lattru2 = 91. ! latsw = gfld%igdtmpl(10) ! lonsw = gfld%igdtmpl(11) ELSE WRITE(6,'(1x,a,I3)') 'Unkown projection: ',gfld%igdtnum CALL arpsstop('Unknown map project in GRIB2 file.',1) END IF CALL gf_free( gfld ) !----------------------------------------------------------------------- ! ! Close file and release the gribfield before returning ! !----------------------------------------------------------------------- CALL baclose( grbunit,iret ) CALL retunit( grbunit ) DEALLOCATE (cgrib) RETURN END SUBROUTINE rdgrb2dims