!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRIBDEC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gribdec(npts,nlev, & 3,12
wdlen,ipds,igds,ibdshd,nbufsz,mgrib, &
fvar,ivar, bds,ibds)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make a GRIB message for variable fvar or ivar for given PDS and
! GDS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! npts Number of grid points
! nlev Number of levels
!
! wdlen Length of machine word
!
! ipds Integer array of PDS
! igds Integer array of GDS
! ibdshd Integer array of BDS header
!
! nbufsz Buffer size of GRIB message
! msglen Length of the GRIB message
! mgrib Buffer carrying the GRIB message
!
! OUTPUT:
!
! pds PDS (GRIB Section 1).
! gds GDS (GRIB Section 3).
! bds BDS except the first 11 octets header
! ibds BDS's integer array
!
! fvar Floating array to be written into GRIB file
! ivar Integer array to be written into GRIB file
!
! TEMPORARY:
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: npts ! Number of grid points
INTEGER :: nlev ! Number of levels
INTEGER :: wdlen ! Length of machine word in bytes
INTEGER :: ipds(25)
INTEGER :: igds(25)
INTEGER :: ibdshd(4)
INTEGER :: nbufsz ! Size of GRIB message array
CHARACTER (LEN=1) :: mgrib(nbufsz) ! GRIB message array
REAL :: fvar(npts) ! Floating array
INTEGER :: ivar(npts) ! Integer array
INTEGER :: ibds(nbufsz/4) ! Identical to BDS
CHARACTER (LEN=1) :: bds(nbufsz) ! BDS
!
!-----------------------------------------------------------------------
!
! Local GRIB parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: bdslen
INTEGER :: pdslen
INTEGER :: gdslen
INTEGER :: msglen ! Length of GRIB message
CHARACTER (LEN=1) :: pds(28) ! PDS
CHARACTER (LEN=1) :: gds(42) ! GDS
CHARACTER (LEN=1) :: bdshd(11) ! BDS header
INTEGER :: ipdsin(25)
INTEGER :: igdsin(25)
INTEGER :: ibdshdin(4)
INTEGER :: iver,bwidth,bscl
INTEGER :: itema
REAL :: vref, scale
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i
INTEGER :: npos, iskip1, iskip2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
npos = 0
iskip1 = 0
iskip2 = 0
!
!-----------------------------------------------------------------------
!
! Check if the record is a GRIB message
!
!-----------------------------------------------------------------------
!
iver = ICHAR(mgrib(8))
IF ( mgrib(1) /= 'G' .OR. mgrib(2) /= 'R' .OR. &
mgrib(3) /= 'I' .OR. mgrib(4) /= 'B' ) THEN
WRITE (6,'(a/a)') &
'Error: this record is not GRIB message', &
'Program stopped in GRIBDEC'
CALL arpsstop
('arpsstop called from gribdec mgrib error',1)
ELSE IF ( iver /= 1 ) THEN
WRITE (6,'(a,i1/a)') &
'Error: GRIB version ID was not correct, iver = ',iver, &
'Program stopped in GRIBDEC'
CALL arpsstop('arpsstop called from gribdec iver error ',1)
END IF
msglen = ICHAR(mgrib(5))*65536 &
+ ICHAR(mgrib(6))*256 &
+ ICHAR(mgrib(7))
npos = npos + 8
!
!-----------------------------------------------------------------------
!
! Get IPDS - GRIB Section 1
!
!-----------------------------------------------------------------------
!
pdslen = ICHAR(mgrib(npos+1))*65536 &
+ ICHAR(mgrib(npos+2))*256 &
+ ICHAR(mgrib(npos+3))
IF ( pdslen /= 28 ) THEN
WRITE (6,'(a,i5/a)') &
'PDS length is not equal to 28, pdslen = ',pdslen, &
'Program stopped in GRIBDEC.'
CALL arpsstop('arpsstop called from gribdec pdslen error',1)
END IF
DO i=1,pdslen
pds(i) = mgrib(npos+i)
END DO
npos = npos + pdslen
CALL gtipds( pds, ipdsin )
IF ( ipdsin(8) /= ipds(8) .OR. ipdsin(9) /= ipds(9) ) THEN
WRITE (6,'(a)') 'ERROR, variable mismatch (ipds 8 or 9).'
GO TO 150
END IF
IF ( ipds(25) /= 0 ) THEN
WRITE (6,'(a,i6,a)') &
'Incorrect decimal scale factor, D-scale = ',ipds(25), &
'ARPS GRIB dump should always has zero D-scale.'
GO TO 150
END IF
DO i=1,25
IF ( ipdsin(i) /= ipds(i) ) THEN
WRITE (6,'(a,i2)') &
'WARNING: A difference was found in IPDS(i): i = ',i
END IF
END DO
! 120 GO TO 170
GO TO 170
150 CONTINUE
WRITE (6,'(a)') &
' i ipds ipdsin'
DO i=1,25
WRITE (6,'(i2,4x,i8,4x,i8)') i, ipds(i),ipdsin(i)
END DO
WRITE (6,*) 'ERROR, critical mismatch found in ', &
'grib file, subroutine GRIBDEC, ABORTING.'
CALL arpsstop('arpsstop called from gribdec file mismatch error',1)
170 CONTINUE
!
!-----------------------------------------------------------------------
!
! Get IGDS - GRIB Section 2
!
!-----------------------------------------------------------------------
!
gdslen = ICHAR(mgrib(npos+1))*65536 &
+ ICHAR(mgrib(npos+2))*256 &
+ ICHAR(mgrib(npos+3))
DO i=1,gdslen
gds(i) = mgrib(npos+i)
END DO
npos = npos + gdslen
CALL gtigds( gds,igdsin )
! DO 210 i=1,5
! IF ( igdsin(i).ne.igds(i) ) THEN
! write (6,'(a,i2)')
! : 'Error: A difference was found in IGDS(i): i = ',i
! GOTO 230
! ENDIF
!210 CONTINUE
GO TO 250
! 230 CONTINUE
! DO i=1,25
! WRITE (6,'(i2,4x,i8,4x,i8)') i, igds(i),igdsin(i)
! END DO
! WRITE (6,'(a)') 'Program stopped in GRIBDEC.'
! CALL arpsstop("arpsstop called from gribdec at 230",1)
250 CONTINUE
! DO 220 i=6,25
! IF ( igdsin(i).ne.igds(i) ) THEN
! write (6,'(a/a)')
! : 'Warning: Data read in may contain errors.',
! : 'A difference was found between IGDS and read-in IGDS:'
! write (6,'(i2,4x,i8,4x,i8)') i, igds(i),igdsin(i)
! ENDIF
!220 CONTINUE
!
!-----------------------------------------------------------------------
!
! Making BDS - GRIB Section 4
!
!-----------------------------------------------------------------------
!
DO i=1,11
bdshd(i) = mgrib(npos+i)
END DO
npos = npos + 11
bdslen = ICHAR(bdshd(1))*65536 &
+ ICHAR(bdshd(2))*256 &
+ ICHAR(bdshd(3))
itema = ICHAR(bdshd(4))
ibdshdin(1) = ishft(IAND(itema,128),-7)
ibdshdin(2) = ishft(IAND(itema, 64),-6)
ibdshdin(3) = ishft(IAND(itema, 32),-5)
ibdshdin(4) = ishft(IAND(itema, 16),-4)
DO i=1,4
IF ( ibdshdin(i) /= ibdshd(i) ) THEN
WRITE (6,'(a,i2.2,a,i4,a,i2.2,a,i4/a)') &
'Error in read-in IBDSHD: ibdshdin(',i,') = ', ibdshdin(i), &
', ibdshd(',i,') = ', ibdshd(i), &
'Program stopped in GRIBDEC'
CALL arpsstop('arpsstop called from gribdec in reading ibdshd',1)
END IF
END DO
bscl = ICHAR(bdshd(5))*256 + ICHAR(bdshd(6))
IF ( bscl > 32768 ) bscl = 32768 - bscl
IF ( ibdshd(3) == 1 .AND. bscl /= 0 ) THEN
WRITE (6,'(a,a,i2,a,i6)') &
'The binary scale factor should be zero for ', &
'original integer data. itype = ',ibdshd(3),', b-scale = ',bscl, &
'Program stopped in GRIBDEC'
CALL arpsstop('arpsstop called from gribdec with bscl ',1)
END IF
scale = 2.**bscl
CALL ibm2flt( bdshd(7), vref )
! write (6,'(a,e20.12)')
! :'The reference value is Vref = ',vref
bwidth = ICHAR(bdshd(11))
DO i=1,bdslen-11
bds(i) = mgrib(npos+i)
END DO
npos = npos + bdslen - 11
IF ( bwidth /= 0 ) THEN
CALL grbgbytes( ibds,ivar,iskip1,bwidth,iskip2, npts )
ELSE
DO i=1,npts
ivar(i) = 0
END DO
END IF
IF ( ibdshd(3) == 0 ) THEN
DO i=1,npts
fvar(i) = FLOAT(ivar(i)) * scale + vref
END DO
ELSE
DO i=1,npts
ivar(i) = ivar(i) + nint(vref)
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Check the end of GRIB message: '7777'
!
!-----------------------------------------------------------------------
!
DO i=1,4
IF ( mgrib(npos+i) /= '7' ) THEN
WRITE (6,'(a)') &
'Incorrect GRIB end flag: ',mgrib(npos+i), &
'It should be alway character ''7''.', &
'Program stopped in GRIBDEC'
CALL arpsstop('arpsstop called from gribdec message 7777',1)
END IF
END DO
npos = npos + 4
IF ( npos /= msglen ) THEN
WRITE (6,'(a/a,i8,a,i8)') &
'The GRIB message length was incorrect.', &
'npos = ',npos, ' msglen = ',msglen
CALL arpsstop('arpsstop called from gribdec with msglen',1)
END IF
RETURN
END SUBROUTINE gribdec
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDHISHD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdhishd(nx,ny,nz,inch,grdbas, fmtver, time, x,y,z, & 1,20
nbufsz,mgrib)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read the header of ARPS history dump by using GRIB representation
! to floating point numbers.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/29/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! inch I/O unit of the history file
!
! grdbas Flag for grid and base state dump
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
!
! fmtver Format and version identical
!
! nbufsz Buffer size of the GRIB message
!
! OUTPUT:
!
! mgrib String containing the header of ARPS history dump
!
! WORK ARRAY:
!
!-----------------------------------------------------------------------
!
! The following parameters are passed into this subroutine through
! a common block in globcst.inc, and they determine which
! variables are output.
!
! grdout =0 or 1. If grdout=0, grid variables are not dumped.
! basout =0 or 1. If basout=0, base state variables are not dumped.
! varout =0 or 1. If varout=0, model perturbation variables are not dumped.
! mstout =0 or 1. If mstout=0, water variables are not dumped.
! rainout=0 or 1. If rainout=0, rain variables are not dumped.
! prcout =0 or 1. If prcout=0, precipitation rates are not dumped.
! iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped.
! tkeout =0 or 1. If tkeout=0, tke is not dumped.
! trbout =0 or 1. If trbout=0, kmh and kmv are not dumped.
! sfcout =0 or 1. If sfcout=0, surface variables are not dumped.
! landout=0 or 1. If landout=0, surface propertty arrays are not dumped.
! ...
!
! These following parameters are also passed in through common
! blocks in globcst.inc.
!
! runname,curtim,umove,vmove,xgrdorg,ygrdorg
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: inch ! I/O unit
INTEGER :: grdbas ! Flag for grid and base state dump
CHARACTER (LEN=40) :: fmtverin
REAL :: x(nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y(ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z(nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
CHARACTER (LEN=40) :: fmtver
INTEGER :: nbufsz
CHARACTER (LEN=1) :: mgrib(nbufsz) ! Header of ARPS GRIB file
REAL :: time
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nxin,nyin,nzin
INTEGER :: npos,hhdlen
INTEGER :: i,j,k
! integer idummy
! real rdummy
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'indtflg.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ (inch) (mgrib(i),i=1,3)
npos = 0
hhdlen = ICHAR(mgrib(npos+1))*65536 &
+ ICHAR(mgrib(npos+2))*256 &
+ ICHAR(mgrib(npos+3))
BACKSPACE (inch)
READ (inch) (mgrib(i),i=1,hhdlen)
npos = npos + 3
DO i=1,40
fmtverin(i:i) = mgrib(npos+i)
END DO
IF ( fmtverin /= fmtver ) THEN
WRITE (6,'(/1x,a,/1x,2a,/1x,3a)') &
'Data format incompatible with the data reader.', &
'Format of data is ',fmtverin,' Format of reader is ',fmtver, &
'. Job stopped.'
CALL arpsstop('arpsstop called from rdhishd with fmtver',1)
END IF
npos = npos + 40
DO i=1,80
runname(i:i) = mgrib(npos+i)
END DO
npos = npos + 80
nocmnt = ICHAR(mgrib(npos+1))
npos = npos + 1
IF( nocmnt > 0 ) THEN
DO j=1,nocmnt
DO i=1,80
cmnt(j)(i:i) = mgrib(npos+i)
END DO
npos = npos + 80
END DO
END IF
CALL ibm2flt( mgrib(npos+1), time )
npos = npos + 4
nxin = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
nyin = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
nzin = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE (6,'(1x,a/1x,3(a,I4)/1x,3(a,I4)/1x,a)') &
'Dimensions in RDHISHD inconsistent with data.', &
'nx = ', nx, ' ny = ', ny, ' nz = ', nz, &
'nxin = ', nxin, ' nyin = ', nyin, ' nzin = ', nzin, &
'Program aborted in RDHISHD.'
CALL arpsstop('arpsstop called from rdhishd with nxin...',1)
END IF
grdin = ICHAR(mgrib(npos+1))
basin = ICHAR(mgrib(npos+2))
varin = ICHAR(mgrib(npos+3))
mstin = ICHAR(mgrib(npos+4))
icein = ICHAR(mgrib(npos+5))
trbin = ICHAR(mgrib(npos+6))
sfcin = ICHAR(mgrib(npos+7))
rainin = ICHAR(mgrib(npos+8))
landin = ICHAR(mgrib(npos+9))
totin = ICHAR(mgrib(npos+10))
tkein = ICHAR(mgrib(npos+11))
mapproj = ICHAR(mgrib(npos+14))
month = ICHAR(mgrib(npos+15))
day = ICHAR(mgrib(npos+16))
year = ICHAR(mgrib(npos+17))*256 + ICHAR(mgrib(npos+18))
hour = ICHAR(mgrib(npos+19))
minute = ICHAR(mgrib(npos+20))
second = ICHAR(mgrib(npos+21))
npos = npos + 21
CALL ibm2flt( mgrib(npos+1), umove )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), vmove )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), xgrdorg )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), ygrdorg )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), trulat1 )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), trulat2 )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), trulon )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), sclfct )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
! CALL ibm2flt( mgrib(npos+1), rdummy )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), tstop )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), thisdmp )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), latitud )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), ctrlat )
npos = npos + 4
CALL ibm2flt( mgrib(npos+1), ctrlon )
npos = npos + 4
IF ( totin /= 0 ) THEN
nstyp = ICHAR(mgrib(npos+1))
IF ( nstyp < 0 ) THEN
nstyp = 1
END IF
npos = npos + 1
!
!-----------------------------------------------------------------------
!
! Reserve 20 integers for future compatibility.
!
!-----------------------------------------------------------------------
!
prcin = ICHAR(mgrib(npos+1))
npos = npos + 1
radin = ICHAR(mgrib(npos+1))
npos = npos + 1
flxin = ICHAR(mgrib(npos+1))
npos = npos + 1
snowcin = ICHAR(mgrib(npos+1))
npos = npos + 1
snowin = ICHAR(mgrib(npos+1))
npos = npos + 1
DO i=6,20
! idummy = ichar(mgrib(npos+1))
npos = npos + 1
END DO
DO i=0,16,4
! CALL ibm2flt( mgrib(npos+i+1), rdummy )
npos = npos + 4
END DO
END IF
!
!-----------------------------------------------------------------------
!
! If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
IF( grdin == 1 .OR. grdbas == 1 ) THEN
DO i=1,nx
CALL ibm2flt( mgrib(npos+1), x(i) )
npos = npos + 4
END DO
DO j=1,ny
CALL ibm2flt( mgrib(npos+1), y(j) )
npos = npos + 4
END DO
DO k=1,nz
CALL ibm2flt( mgrib(npos+1), z(k) )
npos = npos + 4
END DO
END IF
IF ( npos /= hhdlen ) THEN
WRITE (6,'(1x,a,i6,a,i6/a)') &
'Length of header was incorrect. npos = ',npos, &
', hhdlen = ',hhdlen, &
'Program stopped in GTHISHD'
CALL arpsstop('arpsstop called from ghishd with npos ',1)
END IF
RETURN
END SUBROUTINE rdhishd
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE IBM2FLT ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE ibm2flt( cibm, flt ) 19
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Convert GRIB IBM floating representation to machine floating
! number.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/29/1995
! MODIFICATIONS:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! cibm 4 bytes character string containing GRIB IBM format
! floating representation
!
! OUTPUT:
!
! flt Machine dependent floating number
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=1) :: cibm(4)
REAL :: flt
!
!-----------------------------------------------------------------------
!
! Local variables
!
!-----------------------------------------------------------------------
!
REAL :: tema
INTEGER :: ISIGN
INTEGER :: iexp,imant
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
iexp = ICHAR(cibm(1))
IF ( iexp > 127 ) THEN
iexp = iexp - 128
ISIGN = -1
ELSE
ISIGN = 1
END IF
imant = ICHAR(cibm(2))*65536 &
+ ICHAR(cibm(3))*256 &
+ ICHAR(cibm(4))
tema = 16.0**(iexp-70)
flt = FLOAT(ISIGN*imant)*tema
RETURN
END SUBROUTINE ibm2flt
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GTIPDS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gtipds( pds,ipds ) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make PDS from ipds
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/30/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! pds PDS (GRIB Section 1).
!
! OUTPUT:
!
! ipds Integer array of PDS
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=1) :: pds(28) ! PDS
INTEGER :: ipds(25)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Initialize array ipds
!
!-----------------------------------------------------------------------
!
DO i=1,25
ipds(i) = 0
END DO
!
!-----------------------------------------------------------------------
!
! Put pds into ipds
!
!-----------------------------------------------------------------------
!
ipds(1) = ICHAR(pds(1))*65536 &
+ ICHAR(pds(2))*256 &
+ ICHAR(pds(3))
ipds(2) = ICHAR(pds(4))
ipds(3) = ICHAR(pds(5))
ipds(4) = ICHAR(pds(6))
ipds(5) = ICHAR(pds(7))
i = ICHAR(pds(8))
ipds(6) = ishft(IAND(i,128),-7)
ipds(7) = ishft(IAND(i, 64),-6)
ipds(8) = ICHAR(pds(9))
ipds(9) = ICHAR(pds(10))
!
!-----------------------------------------------------------------------
!
! check if level is in two octets or one
!
!-----------------------------------------------------------------------
!
i = ipds(9)
IF ( (i >= 1.AND.i <= 100).OR.i == 102.OR. &
i == 103.OR.i == 105.OR.i == 107.OR. &
i == 109.OR.i == 111.OR.i == 113.OR. &
i == 115.OR.i == 125.OR.i == 160.OR. &
i == 200.OR.i == 201 ) THEN
ipds(11) = ICHAR(pds(11))*256 &
+ ICHAR(pds(12)) ! one level stored two bytes
IF ( ipds(11) > 32768 ) THEN
ipds(11) = 32768 - ipds(11) ! positive integer for the abs(level)
END IF
ELSE
ipds(10) = ICHAR(pds(11))
ipds(11) = ICHAR(pds(12))
END IF
ipds(12) = ICHAR(pds(13))
ipds(13) = ICHAR(pds(14))
ipds(14) = ICHAR(pds(15))
ipds(15) = ICHAR(pds(16))
ipds(16) = ICHAR(pds(17))
ipds(17) = ICHAR(pds(18))
!
!-----------------------------------------------------------------------
!
! Check if time P1 is in two octets or one
!
!-----------------------------------------------------------------------
!
ipds(20) = ICHAR(pds(21))
IF (ipds(20) == 10) THEN
ipds(18) = ICHAR(pds(19))*256 &
+ ICHAR(pds(20))
ELSE
ipds(18) = ICHAR(pds(19))
ipds(19) = ICHAR(pds(20))
END IF
ipds(21) = ICHAR(pds(22))*256 &
+ ICHAR(pds(23))
ipds(22) = ICHAR(pds(24))
ipds(23) = ICHAR(pds(25))
ipds(24) = ICHAR(pds(26))
ipds(25) = ICHAR(pds(27))*256 &
+ ICHAR(pds(28))
IF ( ipds(25) > 32768 ) THEN
ipds(25) = 32768 - ipds(25) ! negative value
END IF
RETURN
END SUBROUTINE gtipds
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GTIGDS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gtigds(gds,igds) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make GDS from igds
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! gds GDS
!
! OUTPUT:
!
! igds Integer array of GDS
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=1) :: gds(42) ! GDS
INTEGER :: igds(25)
!
!-----------------------------------------------------------------------
!
! Local GRIB parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Initialize array igds
!
!-----------------------------------------------------------------------
!
DO i=1,25
igds(i) = 0
END DO
!
!-----------------------------------------------------------------------
!
! Put GDS length in bytes 1,2,3
!
!-----------------------------------------------------------------------
!
igds(1) = ICHAR(gds(4)) ! nz use one byte
igds(2) = ICHAR(gds(5)) ! = 255, N/A for PV & PL
igds(3) = ICHAR(gds(6)) ! data representation
igds(8) = ICHAR(gds(17)) ! resolution and u&v component flag
!
!-----------------------------------------------------------------------
!
! ARPS only support 3 types of projection
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Mercator projection
!
!-----------------------------------------------------------------------
!
IF ( igds(3) == 1 ) THEN
igds(4) = ICHAR(gds(7))*256 + ICHAR(gds(8)) ! nx
igds(5) = ICHAR(gds(9))*256 + ICHAR(gds(10)) ! ny
igds(6) = ICHAR(gds(11))*65536 & ! latitude of 1st point
+ ICHAR(gds(12))*256 &
+ ICHAR(gds(13))
IF ( igds(6) > 8388608 ) THEN
igds(6) = 8388608 - igds(6) ! negative value
END IF
igds(7) = ICHAR(gds(14))*65536 & ! longitude of 1st point
+ ICHAR(gds(15))*256 &
+ ICHAR(gds(16))
IF ( igds(7) > 8388608 ) THEN
igds(7) = 8388608 - igds(7) ! negative value
END IF
igds(9) = ICHAR(gds(18))*65536 & ! latitude of last point
+ ICHAR(gds(19))*256 &
+ ICHAR(gds(20))
IF ( igds(9) > 8388608 ) THEN
igds(9) = 8388608 - igds(9) ! negative value
END IF
igds(10) = ICHAR(gds(21))*65536 & ! longitude of last point
+ ICHAR(gds(22))*256 &
+ ICHAR(gds(23))
IF ( igds(10) > 8388608 ) THEN
igds(10) = 8388608 - igds(10) ! negative value
END IF
igds(11) = ICHAR(gds(32))*65536 & ! d(lat) in milidegree
+ ICHAR(gds(33))*256 &
+ ICHAR(gds(34))
igds(12) = ICHAR(gds(29))*65536 & ! d(lon) in milidegree
+ ICHAR(gds(30))*256 &
+ ICHAR(gds(31))
igds(13) = ICHAR(gds(24))*65536 & ! true latitude
+ ICHAR(gds(25))*256 &
+ ICHAR(gds(26))
igds(14) = ICHAR(gds(28))
!
!-----------------------------------------------------------------------
!
! Lambert Conformal projection
!
!-----------------------------------------------------------------------
!
ELSE IF ( igds(3) == 3 ) THEN
igds(4) = ICHAR(gds(7))*256 & ! nx
+ ICHAR(gds(8))
igds(5) = ICHAR(gds(9))*256 & ! ny
+ ICHAR(gds(10))
igds(6) = ICHAR(gds(11))*65536 &
+ ICHAR(gds(12))*256 &
+ ICHAR(gds(13))
IF ( igds(6) > 8388608 ) THEN
igds(6) = 8388608 - igds(6) ! negative value
END IF
igds(7) = ICHAR(gds(14))*65536 &
+ ICHAR(gds(15))*256 &
+ ICHAR(gds(16))
IF ( igds(7) > 8388608 ) THEN
igds(7) = 8388608 - igds(7) ! negative value
END IF
igds(9) = ICHAR(gds(18))*65536 &
+ ICHAR(gds(19))*256 &
+ ICHAR(gds(20))
IF ( igds(9) > 8388608 ) THEN
igds(9) = 8388608 - igds(9) ! negative value
END IF
igds(10) = ICHAR(gds(21))*65536 & ! dx
+ ICHAR(gds(22))*256 &
+ ICHAR(gds(23))
igds(11) = ICHAR(gds(24))*65536 & ! dy
+ ICHAR(gds(25))*256 &
+ ICHAR(gds(26))
igds(12) = ICHAR(gds(27))
igds(13) = ICHAR(gds(28))
igds(15) = ICHAR(gds(29))*65536 & ! 1st true latitude
+ ICHAR(gds(30))*256 &
+ ICHAR(gds(31))
igds(16) = ICHAR(gds(32))*65536 & ! 2nd true latitude
+ ICHAR(gds(33))*256 &
+ ICHAR(gds(34))
igds(17) = ICHAR(gds(35))*65536 & ! latitude of south pole
+ ICHAR(gds(36))*256 &
+ ICHAR(gds(37))
IF ( igds(17) > 8388608 ) THEN
igds(17) = 8388608 - igds(17) ! negative value
END IF
igds(18) = ICHAR(gds(38))*65536 & ! longitude of south pole
+ ICHAR(gds(39))*256 &
+ ICHAR(gds(40))
IF ( igds(18) > 8388608 ) THEN
igds(18) = 8388608 - igds(18) ! negative value
END IF
!
!-----------------------------------------------------------------------
!
! Polar stereographic projection
!
!-----------------------------------------------------------------------
!
ELSE IF ( igds(3) == 5 ) THEN
igds(4) = ICHAR(gds(7))*256 & ! nx
+ ICHAR(gds(8))
igds(5) = ICHAR(gds(9))*256 & ! ny
+ ICHAR(gds(10))
igds(6) = ICHAR(gds(11))*65536 & ! latitude of 1st point
+ ICHAR(gds(12))*256 &
+ ICHAR(gds(13))
IF ( igds(6) > 8388608 ) THEN
igds(6) = 8388608 - igds(6) ! negative value
END IF
igds(7) = ICHAR(gds(14))*65536 & ! longitude of 1st point
+ ICHAR(gds(15))*256 &
+ ICHAR(gds(16))
IF ( igds(7) > 8388608 ) THEN
igds(7) = 8388608 - igds(7) ! negative value
END IF
igds(9) = ICHAR(gds(18))*65536 & ! true longitude
+ ICHAR(gds(19))*256 &
+ ICHAR(gds(20))
IF ( igds(9) > 8388608 ) THEN
igds(9) = 8388608 - igds(9) ! negative value
END IF
igds(10) = ICHAR(gds(21))*65536 & ! dx
+ ICHAR(gds(22))*256 &
+ ICHAR(gds(23))
igds(11) = ICHAR(gds(24))*65536 & ! dy
+ ICHAR(gds(25))*256 &
+ ICHAR(gds(26))
igds(12) = ICHAR(gds(27))
igds(13) = ICHAR(gds(28))
END IF
RETURN
END SUBROUTINE gtigds