Lambert Conformal Projection and grib2arl

[jhamilton]

Hi,

I am trying to run hysplit with data from the "harmonie" model .. the model runs on
a lambert conformal projection with [in my case] a grid spacing of 2.5km. It is a
limited area model used by the Hirlam/Aladin community.

When I run grib2arl it gives me the following values ..

"...
Input data array from North to South
Input data corner and spacing: 46.83 345.39 -0.00 2.50
Input data array allocation (nlon,nlat): 529 489
..."

When I run chk_file on the output I get ..

"...
Pole pnt lat/lon : 46.8300 225.3900
Reference lat/lon: 0.0000 2.5000
Grid size (km) : 0.0000
Orientation : 0.0000
Tang Lat / Cone : 0.0000
Synch pnt x,y : 1.0000 1.0000
Synch pnt lat/lon: 46.8300 345.3900
Global meteorology grid
..."

This is wrong so I tried setting DLAT/DLON to 2.5/2.5 within the grib2arl code
and then I get

"...
Input data array from North to South
Input data corner and spacing: 46.83 345.39 -2.50 2.50
Input data array allocation (nlon,nlat): 529 489
..."

chk_file now gives ...

"...
Pole pnt lat/lon : 46.8300 225.3900
Reference lat/lon: 2.5000 2.5000
Grid size (km) : 277.9961
Orientation : 0.0000
Tang Lat / Cone : 0.0000
Synch pnt x,y : 1.0000 1.0000
Synch pnt lat/lon:-1000.0000 345.3900
Global meteorology grid
..."

I suspect "harmonie" uses a different lambert grid to the one expected by grib2arl.
The grid is defined by a central lat/long location, a grid apacing [in Km] and a
rotation .. The following is from part of the "harmonie" code ..

"...
C xlov : Longitude of oriention meridian,
C (y increases along this longitude)
C ylat1 : Latitude of first cone intersection latitude
C ylat2 : Latitude of second cone intersection latitude
C xlong : Longitude of origin point of grid
C ylat : Latitude of origin point of grid
C xgrid : Grid spacing in x [metres]
C ygrid : Grid spacing in y [metres]
..."

I understand there are various types of lambert projections and I wonder
if "harmonie" is using a version that is not supported by hysplit.

Thanks

James Hamilton
Met Eireann -- Irish Meteorological Service

[jhamilton]

I am trying to run hysplit using data from the harmonie model. This is
a model developed by the hirlam/aladin community. The version I am
using runs on a 2.5 Km grid centerred at arrrox. 54N 6W and
covering an area of about 1000Km x 1000Km. It uses the lambert
conformal projection.

The grid is defined as follows ..

centre = 233;
dataDate = 20130225;
dataTime = 1200;
stepRange = 24;
startStep = 24;
endStep = 24;
Nx = 529;
Ny = 489;
latitudeOfFirstGridPointInDegrees = 46.834;
longitudeOfFirstGridPointInDegrees = -14.609;
LoVInDegrees = 5;
DxInMetres = 2500;
DyInMetres = 2500;
Latin1InDegrees = 53.5;
Latin2InDegrees = 53.5;

(a) I have tried running grib2arl as follows .. [where
Data_3D_harmonie is a file on harmonie data on pressure
levels and Data_orog_harmonie is the harmonie orography]
This produces a file but hysplit cannot plot it even
when I set the central point to 54N 6W [i.e. 54,-6]. The
plotting program just seems to 'hang'.

When I run grib2arl I get the following MESSAGE file

Command Line Arguments:
Cntr - 60 -80
Grid - 3
Size - 100 100 16
SfcP - 1
Rain - 6
Init - 1
---------------------------------------------------------
Opened file: Data_3D_harmonie
-----------------------------------------------------
SUBROUTINE analyze ...
Input buffer allocation: 517984
First time period found: 13 2 25 12
GRIB code tables set for data center: 233 HAR
Vertical levels in kgds(19): 132
Vertical levels from levels: 66
Vertical level allocation: 132
New level found: 1 1000.0000
Upper level variable match: HGTS
Input data array from North to South
Input data corner and spacing: 46.83 345.39 -0.00 2.50
Input data array allocation (nlon,nlat): 529 489
Bit map array allocation: 258681
Upper level variable match: TEMP
Upper level variable match: UWND
Upper level variable match: VWND
Upper level variable match: SPHU
New level found: 2 925.00000
etc......

Running chk_file on DATA.ARL file gives ..

File start time: 13 2 25 12 0
File ending time: 13 2 26 18 0
Last record number : 868
Record length bytes: 258731
Meteo data model : HARG
Grid size x,y,z : 529 489 6
Vertical coordinate: 2
First forecast hour: 0
Last forecast hour : 30
Records per time : 28
Minutes between : 60
Pole pnt lat/lon : 46.8300 225.3900
Reference lat/lon: 0.0000 2.5000
Grid size (km) : 0.0000
Orientation : 0.0000
Tang Lat / Cone : 0.0000
Synch pnt x,y : 1.0000 1.0000
Synch pnt lat/lon: 46.8300 345.3900
Global meteorology grid
6 500.0000 5 HGTS 51 TEMP 201 UWND 127 VWND 124 SPHU 24
5 700.0000 5 HGTS 27 TEMP 136 UWND 67 VWND 161 SPHU 211
4 850.0000 5 HGTS 224 TEMP 16 UWND 15 VWND 217 SPHU 89
3 925.0000 5 HGTS 78 TEMP 234 UWND 106 VWND 180 SPHU 28
2 1000.0000 5 HGTS 1 TEMP 236 UWND 225 VWND 92 SPHU 183
1 0.0000 2 SHGT 0 SHGT 59
etc ..

(b) I next tried
grib2arl -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
and I got an output file identical to the file produced in (a)

(c) Then I tried
grib2arl -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
and again I got identical results

(d) Next I looked at chaging the value of -g in the argument list
to grib2arl. Withe the -g0 setting I got ..
grib2arl -g0 -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
ERROR: ARL grid beyond the input data domain
Output grid: 1 1 2.8181751 320.36975
Input grid: 134.99150 +Infinity 529 489
STOP 104

(e) With the -g1 setting
grib2arl -g1 -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
ERROR: ARL grid beyond the input data domain
Output grid: 1 1 -20.825668 235.00000
Input grid: 100.84361 +Infinity 529 489
STOP 104

(f) With the -g2 setting
grib2arl -g2 -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
ERROR: ARL grid beyond the input data domain
Output grid: 1 1 20.825668 145.00000
Input grid: 64.843605 +Infinity 529 489
STOP 104

(g) With the -g3 setting
grib2arl -g3 -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
This works but gives the same result as in (a)

(h) With the -g4 setting
grib2arl -g4 -n100:100 -x-6.0 -y52.0 -iData_3D_harmonie -cData_orog_harmonie
Segmentation fault (core dumped)

I cannot think of any other settings on grib2arl that I can try. I would
appreciate help on how to get the program running.

Thanks
James Hamilton
Met Eireann -- Irish Meteological Service

[roland.draxler]

The grib2arl program only handles global lat-lon input files. The attached program (amps2arl) is a simplified version that will convert GRIB1 files on a conformal projection to the HYSPLIT format without altering the map projection. However, it is hardwired for a particular application and you will need to customize the grid size and perhaps the variable selection. It needs to be compiled with a link to the hysplit library.

!################################################################################
! convert GRIB1 format data to ARL format withouth changing the projection
! currently handles only lambert and polar stereographic projections
!-------------------------------------------------------------------------------
! Last Revised: 22 Jan 2013
!################################################################################

PROGRAM AMPS2ARL

LOGICAL ftest
CHARACTER FNAME*80
INTEGER*4 handle, fcopen

! check for command line arguments
narg=IARGC()
if(narg.eq.0)then
write(*,*)'Usage: amps2arl [file_1] [file_2] [...]'
write(*,*)'Convert GRIB1 to ARL format without changing map projection'
stop
end if

! process each data file
do while (narg.gt.0)
CALL GETARG(narg,fname)
INQUIRE(file=fname,exist=ftest)
if(ftest)then
HANDLE=FCOPEN(fname,'r')
write(*,*)'Started processing: ',fname
call xtract(handle)
CALL FCCLOS(handle,*900)
900 continue
else
write(*,*)'File not found:',fname
end if
narg=narg-1
end do

! close out time period and write index record
call pakndx(20)
CLOSE (20)
END

!======================================================================================

SUBROUTINE xtract(handle)

! number of levels, variables, and array limits
PARAMETER (NLVL=18, MVAR=6, MAXX=666, MAXY=627, MAXB=1000000)
PARAMETER (MAXC=600000)

! arrays to hold grib, character, and level variable information
INTEGER SIGL(NLVL), NVAR(NLVL), HANDLE
INTEGER VGRIB, VGRIB0(MVAR), VGRIB1(MVAR),SIG0(MVAR), STYP(MVAR)
CHARACTER*4 VCHAR, VCHARP, VCHAR0(MVAR), VCHAR1(MVAR)
REAL CNVRT, CNVRT0(MVAR), CNVRT1(MVAR)

! define arrays for product definition section, grid description,
! and other grib information
INTEGER KPDS(25),KGDS(25),KPTR(25)

! input data buffer and bit map section
CHARACTER BUFF(MAXB)*1
LOGICAL*1 KBMS(MAXB)
LOGICAL FTEST

! unpacked output array
REAL RVAR(MAXX,MAXY),PVAR(MAXX,MAXY)

! packed output array
CHARACTER CVAR(MAXC)*1, FNAME*80

! default information (grid id, output records)
DATA IG/99/, KREC/0/

! predefine number of variables surface and aloft
! first dimension is number of surface variables,
! second dimension is number of levels above surface * number of upper fields
DATA NVAR / 6, 17*5 /

! set output structure in ARL character format
DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M'/
DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' '/

! set output structure by defining GRIB variables (sfc-0 + level-1)
DATA VGRIB0/ 1, 7, 11, 52, 33, 34/
DATA VGRIB1/ 7, 11, 33, 34, 52, 999/

! special surface variables level identification
DATA SIG0/ 0, 0, 2, 2, 10, 10/
DATA STYP/ 1, 1, 105, 105, 105, 105/

! units conversion factors before converting to ARL format
DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0/
DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0/

! set requested pressure levels for input and output
DATA SIGL/0,1000,975,950,925,900,875,850,800,700,600,500,400,300,250,200,150,100/

! remap array from one to two dimensions
REAL SVAR(MAXX*MAXY)
INTEGER SHAPE(2)
DATA SHAPE/MAXX,MAXY/

!-------------------------------------------------------------------------------
! When dealing with some F90 compilers, replace ICHAR below with JCHAR function
CHARACTER(1) :: mychr
INTEGER :: jchar
JCHAR(MYCHR)=IAND(ICHAR(MYCHR),255)
!-------------------------------------------------------------------------------

SAVE krec

!-------------------------------------------------------------------------------
INTERFACE
SUBROUTINE W3FI63(BUFF,KPDS,KGDS,KBMS,RVAR,KPTR,KRET)
CHARACTER(1), INTENT(IN) :: BUFF(:)
LOGICAL(1), INTENT(OUT) :: KBMS(:)
INTEGER, INTENT(OUT) :: KPDS(:)
INTEGER, INTENT(OUT) :: KGDS(:)
REAL, INTENT(OUT) :: RVAR(:)
INTEGER, INTENT(OUT) :: KPTR(:)
INTEGER, INTENT(OUT) :: KRET
END SUBROUTINE w3fi63
END INTERFACE
!-------------------------------------------------------------------------------

! input grib record byte counter
kbyte=0

! read the indicator section
100 call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,8,*900)

! determine the length of the entire grib record
klen=jchar(buff(7))+jchar(buff(6))*256+jchar(buff(5))*65536
if(klen.gt.maxb)then
write(*,*)'Grib record: ',klen
write(*,*)'Exceedes buffer: ',maxb
stop
end if

! load the indicator section and pds segment into the buffer
call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,36,*900)

! product definition section (+8 byte indicator section offset)
koff=8
kvarb=jchar(buff(koff+9))
ktype=jchar(buff(koff+10))
level=jchar(buff(koff+12))+jchar(buff(koff+11))*256
! check the hours at octet 19 and 20 for precipitation summing
ktim1=jchar(buff(koff+19))
ktim2=jchar(buff(koff+20))
! check the octect 21 to remove all averaged fields
ktim3=jchar(buff(koff+21))

! check if 2d variable present in selection table
kl=1
do kv=1,nvar(kl)
vgrib=vgrib0(kv)
vchar=vchar0(kv)
cnvrt=cnvrt0(kv)
! matches id and special level indicator
IF(KVARB.EQ.VGRIB.AND.LEVEL.EQ.SIG0(KV).AND.KTYPE.EQ.STYP(KV)) &
THEN
! Test for averaging period (want no averaging)
! Precipitation (61,63) is totaled for 12 hours and has ktim3 ne 0, so don't check it
IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63)GO TO 300
END IF
end do

! then check for 3d variable
kl=2
do kv=1,nvar(kl)
vgrib=vgrib1(kv)
vchar=vchar1(kv)
cnvrt=cnvrt1(kv)
if(kvarb.eq.vgrib)go to 200
end do
go to 800

! check if 3d level is present in selection table
200 do kl=2,nlvl
if(level.eq.sigl(kl))go to 300
end do

! if all tests fail go and read next grib record
go to 800

! load the entire grib data record into the buffer
300 krec=krec+1
call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,klen,*900)

! call the nmc grib unpacker
call W3FI63(buff,kpds,kgds,kbms,svar,kptr,kret)
! remap input data from one- to two-dimensional array
RVAR=RESHAPE(SVAR,SHAPE)
if(kret.ne.0)then
write(*,*)'Error W3FI63: ',kret
stop
end if

! correct for Y2K
kpds(8)=mod(kpds(8),100)

! set the current time
iyr=kpds(8)
imo=kpds(9)
ida=kpds(10)
ihr=kpds(11)
imn=kpds(12)
ifh=max(kpds(14),kpds(15))
call tmplus(iyr,imo,ida,ihr,ifh)

! set previous time the first time
IF(krec.EQ.1)THEN
kyr=iyr
kmo=imo
kda=ida
khr=ihr

! check for time change
ELSE
IF(iyr.NE.kyr.OR.imo.NE.kmo.OR.ida.NE.kda.OR.ihr.NE.khr)THEN
CALL PAKNDX(20)
WRITE(*,*)'Finished: ',kyr,kmo,kda,khr
kyr=iyr
kmo=imo
kda=ida
khr=ihr
END IF
END IF

! for the the first record create an index record for pakrec
IF(krec.eq.1)then

! decode grib information and create index record structure
call makndx(kgds,vchar0,vchar1,mvar,nlvl,nvar,sigl,ig,nxp,nyp)

if(nxp.NE.maxx.OR.nyp.NE.maxy)then
write(*,*)'Input data array size: ',nxp,nyp
write(*,*)'Compiled dimensions: ',maxx,maxy
stop
end if

nxy=nxp*nyp
if(nxy.gt.maxc)then
write(*,*)'Packed output array: ',nxy
write(*,*)'Exceeds dimensions: ',maxc
stop
end if
call pakset(20,'CFG.ARL',1,nxp,nyp,nzp)

! write control file with date
OPEN(40,FILE='ARLTIME')
WRITE(40,'(5I2)')IYR,IMO,IDA,IHR,IFH
CLOSE (40)

! standard file name for output
FNAME='DATA.ARL'
INQUIRE(FILE=FNAME,EXIST=FTEST)

! open data file
LREC=NXY+50
OPEN(20,FILE=FNAME,RECL=LREC, &
ACCESS='DIRECT',FORM='UNFORMATTED')

IF(.NOT.FTEST)THEN
! output data set initialized to missing
KINI=1
ELSE
! output data set created previously, fill in data
KINI=0
END IF
write(*,*)'Initialized output data set: ',fname
END IF


! precipitation (A PCP - 61, ACPCP = 63) exist at forecast hour +0,
! but don't want to overwrite the +6 hour forecast from the previous cycle
IF(KINI.EQ.0.AND.(VGRIB.EQ.61.OR.VGRIB.EQ.63))GO TO 800

! perform any required units conversion
if(cnvrt.ne.1.0) call datcnv(rvar,nxp,nyp,cnvrt)

! then pack into ARL format and continue
! write(*,*)'Record: ',krec,kv,kl,ktype,vchar,sigl(kl)
call pakrec(20,rvar,cvar,nxp,nyp,nxy,vchar,IYR,IMO,IDA,IHR,IMN,IFH,KL,KINI)

800 kbyte=kbyte+klen
GO TO 100

900 RETURN
END

!=================================================================================================

SUBROUTINE MAKNDX(kgds,vchar0,vchar1,mvar,nlvl,nvar,sigl,ig,nxp,nyp)

LOGICAL FTEST
INTEGER kgds(25)
REAL grids(12)

! arrays to hold variable selection information
INTEGER sigl(nlvl), nvar(nlvl)
CHARACTER*4 VCHAR0(MVAR), VCHAR1(MVAR)

! the number of grid points
nxp=kgds(2)
nyp=kgds(3)

! grid orientation
grids(6)=0.0

! lat/lon of grid point 1,1
grids(8)=1.0
grids(9)=1.0
grids(10)=kgds(4)/1000.0
grids(11)=kgds(5)/1000.0
grids(12)=0.0

! defines a lambert conformal projection
if(kgds(1).eq.3)then

! set the pole position and reference lat/lon
grids(1)=90.0
grids(2)=0.0
grids(3)=kgds(12)/1000.0

! reference longitude = grid orientation
grids(4)=kgds(7)/1000.0

! delta=x grid size
grids(5)=kgds(8)/1000.0

! tangent latitude
grids(7)=kgds(12)/1000.0

! defines a polar sterographic projection
elseif(kgds(1).eq.5)then
! set hemisphere (north=1, south=-1)
khem=1
if(kgds(10).eq.128) khem=-1

! set the pole position and reference lat/lon
grids(1)=90.0*khem

! pole longtitude (+180 from cut)
grids(2)=kgds(7)/1000.0

! reference lat/lon (at which grid size specified)
grids(3)=60.0*khem

! reference longitude and grid alignment
grids(4)=grids(2)

! delta=x grid size
grids(5)=kgds(8)/1000.0

! tangent latitude
grids(7)=grids(1)

else
write(*,*)'Undefined projection: ',kgds(1)
stop
end if

! write the packer configuration file
OPEN(30,file='CFG.ARL')
WRITE(30,'(20X,A4)')'AMPS'

! grid number 99 and 1 for sigma coordinate system
! grid number 99 and 2 for pressure coordinate system
WRITE(30,'(20X,I4)') IG, 2

WRITE(30,'(20X,F10.4)')(GRIDS(I),I=1,12)
WRITE(30,'(20X,I4)')nxp,nyp,nlvl

! upper level information
DO NL=1,NLVL
sigma=sigl(nl)
if(nl.eq.1)then
WRITE(30,'(20X,F6.1,I3,32(1X,A4))')sigma,nvar(nl),(vchar0(nv),nv=1,nvar(nl))
else
WRITE(30,'(20X,F6.1,I3,32(1X,A4))')sigma,nvar(nl),(vchar1(nv),nv=1,nvar(nl))
end if
END DO
CLOSE (30)
RETURN
END

!=====================================================================================

SUBROUTINE datcnv(rvar,nxp,nyp,cnvrt)

REAL RVAR(NXP,NYP)

DO J=1,NYP
DO I=1,NXP

rvar(i,j)=rvar(i,j)*cnvrt

END DO
END DO
RETURN
END

[jhamilton]

Many thanks, I have got your program running with the harmonie data [on pressure levels] and I am now adding in rainfall. I will post a copy of the final program whan I have finished modifying it.
James Hamilton

[jhamilton]

Here are details of the modifications I made to the AMPS2ARL tprogram to get it to work with
the Harmonie model ..

**************************************************************************************************

diff harm2arl.f amps2arl.f
2d1
< ! Postby roland.draxler » February 28th, 2013, 5:21 pm
17,28c16
< !ORG! PROGRAM AMPS2ARL
< PROGRAM HARM2ARL
< !
< ! This is a modification of the AMPS2ARL program to deal with the Harmonie model
< ! as used at Met Eireann .. where lines have been changed the original version
< ! has been indicated with the comment '!ORG!. The following changes have been
< ! made ..
< ! a) changes to grid dimensions and number of levels
< ! b) changes to parameter lists and to parameter GRIB codes
< ! c) the GRIB code for Harmonie rainfall is '181'
< ! d) the programs checks the length of the rainfall period and sets the
< ! parameter name to TPP1, TPP3 or TPP6 as appropriate
---
> PROGRAM AMPS2ARL
68,69c56
< !ORG! PARAMETER (NLVL=18, MVAR=6, MAXX=666, MAXY=627, MAXB=1000000)
< PARAMETER (NLVL=6, MVAR=7, MAXX=529, MAXY=489, MAXB=1000000)
---
> PARAMETER (NLVL=18, MVAR=6, MAXX=666, MAXY=627, MAXB=1000000)
99,100c86
< !ORG! DATA NVAR / 6, 17*5 /
< DATA NVAR / 7, 5*5 /
---
> DATA NVAR / 6, 17*5 /
103,106c89,90
< !ORG! DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M'/
< !ORG! DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' '/
< DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M','TPPX'/
< DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' ',' '/
---
> DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M'/
> DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' '/
109,112c93,94
< !ORG! DATA VGRIB0/ 1, 7, 11, 52, 33, 34/
< !ORG! DATA VGRIB1/ 7, 11, 33, 34, 52, 999/
< DATA VGRIB0/ 1, 6, 11, 52, 33, 34, 181/ !'PRSS','SHGT','T02M','RH2M','U10M','V10M','TPPX'
< DATA VGRIB1/ 6, 11, 33, 34, 52, 999, 999/ !'HGTS','TEMP','UWND','VWND','RELH',' '
---
> DATA VGRIB0/ 1, 7, 11, 52, 33, 34/
> DATA VGRIB1/ 7, 11, 33, 34, 52, 999/
115,118c97,98
< !ORG! DATA SIG0/ 0, 0, 2, 2, 10, 10/
< !ORG! DATA STYP/ 1, 1, 105, 105, 105, 105/
< DATA SIG0/ 0, 0, 2, 2, 10, 10, 0/
< DATA STYP/ 105, 105, 105, 105, 105, 105, 105/
---
> DATA SIG0/ 0, 0, 2, 2, 10, 10/
> DATA STYP/ 1, 1, 105, 105, 105, 105/
121,124c101,102
< !ORG! DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0/
< !ORG! DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0/
< DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
< DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0/
---
> DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0/
> DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0/
127,128c105
< !ORG! DATA SIGL/0,1000,975,950,925,900,875,850,800,700,600,500,400,300,250,200,150,100/
< DATA SIGL/0,1000,925,850,700,500/
---
> DATA SIGL/0,1000,975,950,925,900,875,850,800,700,600,500,400,300,250,200,150,100/
197,208d173
< ! new block added at Met Eireann ....................... start
< ! decide if we have 1-hour, 3-hour or 6-hour rainfall
< if(VGRIB.EQ.61.OR.VGRIB.EQ.63.or.vgrib.eq.181) then
< if(VCHAR0(KV).eq.'TPPX') then
< if((ktim2-ktim1).eq.1) VCHAR0(KV)='TPP1' !1-hour rainfall
< if((ktim2-ktim1).eq.3) VCHAR0(KV)='TPP3' !3-hour rainfall
< if((ktim2-ktim1).eq.6) VCHAR0(KV)='TPP6' !6-hour rainfall
< write(*,*)'Setting Rainfall : ',VCHAR0(KV)
< vchar=VCHAR0(KV)
< endif
< endif
< ! new block added at Met Eireann ......................... end
211,213c176
< ! Harmonie precip parameter code is 181
< !ORG! IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63)GO TO 300
< IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63.or.vgrib.eq.181)GO TO 300
---
> IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63)GO TO 300
416,417c379
< !ORG! WRITE(30,'(20X,A4)')'AMPS'
< WRITE(30,'(20X,A4)')'HARM'
---
> WRITE(30,'(20X,A4)')'AMPS'


**************************************************************************************************

Here is a listing of the whole [modified] program

**************************************************************************************************

!
! Postby roland.draxler » February 28th, 2013, 5:21 pm
! The grib2arl program only handles global lat-lon input files. The attached
! program (amps2arl) is a simplified version that will convert GRIB1 files on
! a conformal projection to the HYSPLIT format without altering the map
! projection. However, it is hardwired for a particular application and you
! will need to customize the grid size and perhaps the variable selection.
! It needs to be compiled with a link to the hysplit library.
!
!################################################################################
! convert GRIB1 format data to ARL format withouth changing the projection
! currently handles only lambert and polar stereographic projections
!-------------------------------------------------------------------------------
! Last Revised: 22 Jan 2013
!################################################################################

!ORG! PROGRAM AMPS2ARL
PROGRAM HARM2ARL
!
! This is a modification of the AMPS2ARL program to deal with the Harmonie model
! as used at Met Eireann .. where lines have been changed the original version
! has been indicated with the comment '!ORG!. The following changes have been
! made ..
! a) changes to grid dimensions and number of levels
! b) changes to parameter lists and to parameter GRIB codes
! c) the GRIB code for Harmonie rainfall is '181'
! d) the programs checks the length of the rainfall period and sets the
! parameter name to TPP1, TPP3 or TPP6 as appropriate

LOGICAL ftest
CHARACTER FNAME*80
INTEGER*4 handle, fcopen

! check for command line arguments
narg=IARGC()
if(narg.eq.0)then
write(*,*)'Usage: amps2arl [file_1] [file_2] [...]'
write(*,*)'Convert GRIB1 to ARL format without changing map projection'
stop
end if

! process each data file
do while (narg.gt.0)
CALL GETARG(narg,fname)
INQUIRE(file=fname,exist=ftest)
if(ftest)then
HANDLE=FCOPEN(fname,'r')
write(*,*)'Started processing: ',fname
call xtract(handle)
CALL FCCLOS(handle,*900)
900 continue
else
write(*,*)'File not found:',fname
end if
narg=narg-1
end do

! close out time period and write index record
call pakndx(20)
CLOSE (20)
END

!======================================================================================

SUBROUTINE xtract(handle)

! number of levels, variables, and array limits
!ORG! PARAMETER (NLVL=18, MVAR=6, MAXX=666, MAXY=627, MAXB=1000000)
PARAMETER (NLVL=6, MVAR=7, MAXX=529, MAXY=489, MAXB=1000000)
PARAMETER (MAXC=600000)

! arrays to hold grib, character, and level variable information
INTEGER SIGL(NLVL), NVAR(NLVL), HANDLE
INTEGER VGRIB, VGRIB0(MVAR), VGRIB1(MVAR),SIG0(MVAR), STYP(MVAR)
CHARACTER*4 VCHAR, VCHARP, VCHAR0(MVAR), VCHAR1(MVAR)
REAL CNVRT, CNVRT0(MVAR), CNVRT1(MVAR)

! define arrays for product definition section, grid description,
! and other grib information
INTEGER KPDS(25),KGDS(25),KPTR(25)

! input data buffer and bit map section
CHARACTER BUFF(MAXB)*1
LOGICAL*1 KBMS(MAXB)
LOGICAL FTEST

! unpacked output array
REAL RVAR(MAXX,MAXY),PVAR(MAXX,MAXY)

! packed output array
CHARACTER CVAR(MAXC)*1, FNAME*80

! default information (grid id, output records)
DATA IG/99/, KREC/0/

! predefine number of variables surface and aloft
! first dimension is number of surface variables,
! second dimension is number of levels above surface * number of upper fields
!ORG! DATA NVAR / 6, 17*5 /
DATA NVAR / 7, 5*5 /

! set output structure in ARL character format
!ORG! DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M'/
!ORG! DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' '/
DATA VCHAR0/ 'PRSS','SHGT','T02M','RH2M','U10M','V10M','TPPX'/
DATA VCHAR1/ 'HGTS','TEMP','UWND','VWND','RELH',' ',' '/

! set output structure by defining GRIB variables (sfc-0 + level-1)
!ORG! DATA VGRIB0/ 1, 7, 11, 52, 33, 34/
!ORG! DATA VGRIB1/ 7, 11, 33, 34, 52, 999/
DATA VGRIB0/ 1, 6, 11, 52, 33, 34, 181/ !'PRSS','SHGT','T02M','RH2M','U10M','V10M','TPPX'
DATA VGRIB1/ 6, 11, 33, 34, 52, 999, 999/ !'HGTS','TEMP','UWND','VWND','RELH',' '

! special surface variables level identification
!ORG! DATA SIG0/ 0, 0, 2, 2, 10, 10/
!ORG! DATA STYP/ 1, 1, 105, 105, 105, 105/
DATA SIG0/ 0, 0, 2, 2, 10, 10, 0/
DATA STYP/ 105, 105, 105, 105, 105, 105, 105/

! units conversion factors before converting to ARL format
!ORG! DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0/
!ORG! DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0/
DATA CNVRT0/ 0.01, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/
DATA CNVRT1/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0/

! set requested pressure levels for input and output
!ORG! DATA SIGL/0,1000,975,950,925,900,875,850,800,700,600,500,400,300,250,200,150,100/
DATA SIGL/0,1000,925,850,700,500/

! remap array from one to two dimensions
REAL SVAR(MAXX*MAXY)
INTEGER SHAPE(2)
DATA SHAPE/MAXX,MAXY/

!-------------------------------------------------------------------------------
! When dealing with some F90 compilers, replace ICHAR below with JCHAR function
CHARACTER(1) :: mychr
INTEGER :: jchar
JCHAR(MYCHR)=IAND(ICHAR(MYCHR),255)
!-------------------------------------------------------------------------------

SAVE krec

!-------------------------------------------------------------------------------
INTERFACE
SUBROUTINE W3FI63(BUFF,KPDS,KGDS,KBMS,RVAR,KPTR,KRET)
CHARACTER(1), INTENT(IN) :: BUFF(:)
LOGICAL(1), INTENT(OUT) :: KBMS(:)
INTEGER, INTENT(OUT) :: KPDS(:)
INTEGER, INTENT(OUT) :: KGDS(:)
REAL, INTENT(OUT) :: RVAR(:)
INTEGER, INTENT(OUT) :: KPTR(:)
INTEGER, INTENT(OUT) :: KRET
END SUBROUTINE w3fi63
END INTERFACE
!-------------------------------------------------------------------------------

! input grib record byte counter
kbyte=0

! read the indicator section
100 call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,8,*900)

! determine the length of the entire grib record
klen=jchar(buff(7))+jchar(buff(6))*256+jchar(buff(5))*65536
if(klen.gt.maxb)then
write(*,*)'Grib record: ',klen
write(*,*)'Exceedes buffer: ',maxb
stop
end if

! load the indicator section and pds segment into the buffer
call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,36,*900)

! product definition section (+8 byte indicator section offset)
koff=8
kvarb=jchar(buff(koff+9))
ktype=jchar(buff(koff+10))
level=jchar(buff(koff+12))+jchar(buff(koff+11))*256
! check the hours at octet 19 and 20 for precipitation summing
ktim1=jchar(buff(koff+19))
ktim2=jchar(buff(koff+20))
! check the octect 21 to remove all averaged fields
ktim3=jchar(buff(koff+21))

! check if 2d variable present in selection table
kl=1
do kv=1,nvar(kl)
vgrib=vgrib0(kv)
vchar=vchar0(kv)
cnvrt=cnvrt0(kv)
! matches id and special level indicator
IF(KVARB.EQ.VGRIB.AND.LEVEL.EQ.SIG0(KV).AND.KTYPE.EQ.STYP(KV)) &
THEN
! new block added at Met Eireann ....................... start
! decide if we have 1-hour, 3-hour or 6-hour rainfall
if(VGRIB.EQ.61.OR.VGRIB.EQ.63.or.vgrib.eq.181) then
if(VCHAR0(KV).eq.'TPPX') then
if((ktim2-ktim1).eq.1) VCHAR0(KV)='TPP1' !1-hour rainfall
if((ktim2-ktim1).eq.3) VCHAR0(KV)='TPP3' !3-hour rainfall
if((ktim2-ktim1).eq.6) VCHAR0(KV)='TPP6' !6-hour rainfall
write(*,*)'Setting Rainfall : ',VCHAR0(KV)
vchar=VCHAR0(KV)
endif
endif
! new block added at Met Eireann ......................... end
! Test for averaging period (want no averaging)
! Precipitation (61,63) is totaled for 12 hours and has ktim3 ne 0, so don't check it
! Harmonie precip parameter code is 181
!ORG! IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63)GO TO 300
IF(KTIM3.EQ.0.OR.VGRIB.EQ.61.OR.VGRIB.EQ.63.or.vgrib.eq.181)GO TO 300
END IF
end do

! then check for 3d variable
kl=2
do kv=1,nvar(kl)
vgrib=vgrib1(kv)
vchar=vchar1(kv)
cnvrt=cnvrt1(kv)
if(kvarb.eq.vgrib)go to 200
end do
go to 800

! check if 3d level is present in selection table
200 do kl=2,nlvl
if(level.eq.sigl(kl))go to 300
end do

! if all tests fail go and read next grib record
go to 800

! load the entire grib data record into the buffer
300 krec=krec+1
call fcptps(handle,kbyte,*900)
call fcread(handle,buff,1,klen,*900)

! call the nmc grib unpacker
call W3FI63(buff,kpds,kgds,kbms,svar,kptr,kret)
! remap input data from one- to two-dimensional array
RVAR=RESHAPE(SVAR,SHAPE)
if(kret.ne.0)then
write(*,*)'Error W3FI63: ',kret
stop
end if

! correct for Y2K
kpds(8)=mod(kpds(8),100)

! set the current time
iyr=kpds(8)
imo=kpds(9)
ida=kpds(10)
ihr=kpds(11)
imn=kpds(12)
ifh=max(kpds(14),kpds(15))
call tmplus(iyr,imo,ida,ihr,ifh)

! set previous time the first time
IF(krec.EQ.1)THEN
kyr=iyr
kmo=imo
kda=ida
khr=ihr

! check for time change
ELSE
IF(iyr.NE.kyr.OR.imo.NE.kmo.OR.ida.NE.kda.OR.ihr.NE.khr)THEN
CALL PAKNDX(20)
WRITE(*,*)'Finished: ',kyr,kmo,kda,khr
kyr=iyr
kmo=imo
kda=ida
khr=ihr
END IF
END IF

! for the the first record create an index record for pakrec
IF(krec.eq.1)then

! decode grib information and create index record structure
call makndx(kgds,vchar0,vchar1,mvar,nlvl,nvar,sigl,ig,nxp,nyp)

if(nxp.NE.maxx.OR.nyp.NE.maxy)then
write(*,*)'Input data array size: ',nxp,nyp
write(*,*)'Compiled dimensions: ',maxx,maxy
stop
end if

nxy=nxp*nyp
if(nxy.gt.maxc)then
write(*,*)'Packed output array: ',nxy
write(*,*)'Exceeds dimensions: ',maxc
stop
end if
call pakset(20,'CFG.ARL',1,nxp,nyp,nzp)

! write control file with date
OPEN(40,FILE='ARLTIME')
WRITE(40,'(5I2)')IYR,IMO,IDA,IHR,IFH
CLOSE (40)

! standard file name for output
FNAME='DATA.ARL'
INQUIRE(FILE=FNAME,EXIST=FTEST)

! open data file
LREC=NXY+50
OPEN(20,FILE=FNAME,RECL=LREC, &
ACCESS='DIRECT',FORM='UNFORMATTED')

IF(.NOT.FTEST)THEN
! output data set initialized to missing
KINI=1
ELSE
! output data set created previously, fill in data
KINI=0
END IF
write(*,*)'Initialized output data set: ',fname
END IF


! precipitation (A PCP - 61, ACPCP = 63) exist at forecast hour +0,
! but don't want to overwrite the +6 hour forecast from the previous cycle
IF(KINI.EQ.0.AND.(VGRIB.EQ.61.OR.VGRIB.EQ.63))GO TO 800

! perform any required units conversion
if(cnvrt.ne.1.0) call datcnv(rvar,nxp,nyp,cnvrt)

! then pack into ARL format and continue
! write(*,*)'Record: ',krec,kv,kl,ktype,vchar,sigl(kl)
call pakrec(20,rvar,cvar,nxp,nyp,nxy,vchar,IYR,IMO,IDA,IHR,IMN,IFH,KL,KINI)

800 kbyte=kbyte+klen
GO TO 100

900 RETURN
END

!=================================================================================================

SUBROUTINE MAKNDX(kgds,vchar0,vchar1,mvar,nlvl,nvar,sigl,ig,nxp,nyp)

LOGICAL FTEST
INTEGER kgds(25)
REAL grids(12)

! arrays to hold variable selection information
INTEGER sigl(nlvl), nvar(nlvl)
CHARACTER*4 VCHAR0(MVAR), VCHAR1(MVAR)

! the number of grid points
nxp=kgds(2)
nyp=kgds(3)

! grid orientation
grids(6)=0.0

! lat/lon of grid point 1,1
grids(8)=1.0
grids(9)=1.0
grids(10)=kgds(4)/1000.0
grids(11)=kgds(5)/1000.0
grids(12)=0.0

! defines a lambert conformal projection
if(kgds(1).eq.3)then

! set the pole position and reference lat/lon
grids(1)=90.0
grids(2)=0.0
grids(3)=kgds(12)/1000.0

! reference longitude = grid orientation
grids(4)=kgds(7)/1000.0

! delta=x grid size
grids(5)=kgds(8)/1000.0

! tangent latitude
grids(7)=kgds(12)/1000.0

! defines a polar sterographic projection
elseif(kgds(1).eq.5)then
! set hemisphere (north=1, south=-1)
khem=1
if(kgds(10).eq.128) khem=-1

! set the pole position and reference lat/lon
grids(1)=90.0*khem

! pole longtitude (+180 from cut)
grids(2)=kgds(7)/1000.0

! reference lat/lon (at which grid size specified)
grids(3)=60.0*khem

! reference longitude and grid alignment
grids(4)=grids(2)

! delta=x grid size
grids(5)=kgds(8)/1000.0

! tangent latitude
grids(7)=grids(1)

else
write(*,*)'Undefined projection: ',kgds(1)
stop
end if

! write the packer configuration file
OPEN(30,file='CFG.ARL')
!ORG! WRITE(30,'(20X,A4)')'AMPS'
WRITE(30,'(20X,A4)')'HARM'

! grid number 99 and 1 for sigma coordinate system
! grid number 99 and 2 for pressure coordinate system
WRITE(30,'(20X,I4)') IG, 2

WRITE(30,'(20X,F10.4)')(GRIDS(I),I=1,12)
WRITE(30,'(20X,I4)')nxp,nyp,nlvl

! upper level information
DO NL=1,NLVL
sigma=sigl(nl)
if(nl.eq.1)then
WRITE(30,'(20X,F6.1,I3,32(1X,A4))')sigma,nvar(nl),(vchar0(nv),nv=1,nvar(nl))
else
WRITE(30,'(20X,F6.1,I3,32(1X,A4))')sigma,nvar(nl),(vchar1(nv),nv=1,nvar(nl))
end if
END DO
CLOSE (30)
RETURN
END

!=====================================================================================

SUBROUTINE datcnv(rvar,nxp,nyp,cnvrt)

REAL RVAR(NXP,NYP)

DO J=1,NYP
DO I=1,NXP

rvar(i,j)=rvar(i,j)*cnvrt

END DO
END DO
RETURN
END


**************************************************************************************************

Thanks again for all the help
James Hamilton
Met Eireann