program MM5toGRIB 

!
! This utility program is written in free-format Fortran 90.
! It requires a Fortran 90 compiler to compile.
!
! This utility programs converts MM5 format output files
! in GRIB format files
!   
! Source programs are:    - readv3 from MM5 utils
!                         - subroutines emoslib from ECMWF for WMO FM-92 GRIB code
!
! 2004 Instituto Nacional de Meteorologia (Spain)
! Servicio de Predecibilidad y Predicciones Extendidas
! Author: Alfons Callado Pallares
! E-mail: acallado@inm.es
! Last update: 20061214
!

  implicit none

  common /MM5SEC/ ordering, name, staggering, current_date, bhi, time,&
                  PRESSURE, LATITDOT, LATITCRS, LONGIDOT, LONGICRS,&
                  ENDLATITDOT, ENDLATITCRS, ENDLONGIDOT, ENDLONGICRS, bhr,&
                  nfield
  common /interpol/ lweight,lcross,lUV,end_index,end_index_n,&
                    inter,start_lon,start_lat,interval_lon,interval_lat

  integer, dimension(50,20) :: bhi
  real, dimension(20,20) :: bhr
  character(len=80), dimension(50,20) :: bhic
  character(len=80), dimension(20,20) :: bhrc
  character(len=120) :: flnm
  character*40 :: NAMEFILE
  integer*4 :: iunit = 10

  integer :: flag

  integer :: ndim
  real :: time, sample
  integer, dimension(4) :: start_index, end_index,end_index_n
  character (len= 4) :: staggering
  character (len= 4) :: ordering
  character (len=24) :: start_date
  character (len=24) :: current_date
  character (len= 9) :: name
  integer :: ipp
  character (len=25) :: units
  character (len=46) :: description

  integer :: l
  integer :: inn_end
  parameter (inn_end = 200)
  character (len= 9), dimension(inn_end) :: nametogrib
  character (len= 9) :: nameT
  character (len= 9) :: nameRT = 'RAIN TOT '
  logical, dimension(inn_end) :: ltogrib = .false.
  logical :: ltogribyes = .false.
  integer :: inn,innT,innn,innn_end
  integer,dimension(inn_end) :: in

  real, allocatable, dimension(:,:,:,:) :: data
  real, allocatable, dimension(:,:,:,:) :: dataTRain
  real, allocatable, dimension(:,:,:,:) :: datan
  real, allocatable, dimension(:,:,:,:,:) :: datanUV
  real, allocatable, dimension(:,:,:) :: wiscrs
  real, allocatable, dimension(:,:,:) :: wisdot
  real, allocatable, dimension(:,:,:) :: wagdot
  integer, allocatable, dimension(:,:,:,:) :: index
  integer, allocatable, dimension(:,:,:,:) :: indexUV

  real :: start_lon,start_lat,interval_lon,interval_lat
  integer :: inter

  integer :: ierrr, ierr, ier

  logical :: newtime = .TRUE.
  logical :: lmore
  logical :: lsplit
  logical :: loneg = .true.
  logical :: llonlat
  logical :: lweight = .true.
  logical :: lcross
  integer :: lUV = 0

  logical :: lRainTot = .false.
  logical :: lRainNon = .false.
  logical :: lRainCon = .false.
  integer :: lTRain = 0

  integer :: LATITDOT=99999,LATITCRS=99999,LONGIDOT=999999,LONGICRS=999999
  integer :: ENDLATITDOT,ENDLATITCRS,ENDLONGIDOT,ENDLONGICRS
  integer :: NPRESSURE
  parameter (NPRESSURE=50)
  integer, dimension(NPRESSURE) :: PRESSURE=-99999
  real, allocatable, dimension(:,:) :: loncrs, latcrs, londot, latdot 
  logical :: lLALO = .FALSE.

  type sksecone
          integer :: six
          integer :: seven
          integer :: eight
          integer :: nine
          real    :: factor
  end type sksecone

  type(sksecone),dimension(inn_end) :: sksec
  type(sksecone),dimension(inn_end) :: rksec
  type(sksecone),dimension(inn_end) :: rksecU
  type(sksecone),dimension(inn_end) :: rksecV
  type(sksecone),dimension(inn_end) :: rksecRT
  integer :: nfield

  call arguments(flnm, lmore, lsplit, llonlat)

  call read_namelist_output(ltogrib,nametogrib,lmore,sksec)
!    stop

  if (lmore) print*, 'flnm = ', trim(flnm)
  open(iunit, file=flnm, err=99, form='unformatted', status='old', action='read')

  read(iunit, iostat=ierr) flag
  if (lmore) print*, iunit, ierr, flag
  ierrr=ierr
  do while (ierr == 0)

     if (flag == 0) then
        read(iunit,iostat=ier) bhi, bhr, bhic, bhrc
        lweight= .true. 
        if(ier/=0) then
           write(*,'("Error reading big header")')
           call abort()
        endif
        if ((.NOT.(lLALO)).and.(lmore))  then
           call printout_big_header(bhi, bhr, bhic, bhrc)
        endif
     elseif (flag == 1) then

        READ (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering,&
             current_date, name, units, description
        if(ier/=0) then
           write(*,'("Error reading subheader")')
           call abort()
        endif

        if ((lmore).AND.(lLALO)) then
           print*, 'ndim: ', ndim
           print*, 'start_index: ', start_index
           print*, 'end_index: ', end_index
           print*, 'time: ', time
           print*, 'staggering: #'//staggering//'#'
           print*, 'ordering: #'//ordering//'#'
           print*, 'date/time: #'//current_date//'#'
           print*, 'name: #'//name//'#'
           print*, 'units: #'//units//'#'
           print*, 'description: #'//description//'#'
        endif

        if ((newtime).AND.(lLALO)) then
           if (lmore) write(*,'(/,A,2x, F15.5," Hours"/)') current_date, time/60.
           newtime = .FALSE.
        endif

        if (ndim == 1) then
           allocate(data(end_index(1), 1, 1, 1))
        elseif (ndim == 2) then
           allocate(data(end_index(1), end_index(2), 1, 1))
        elseif (ndim == 3) then
           allocate(data(end_index(1), end_index(2), end_index(3), 1))
        endif

        read(iunit) data

!
!       lLALO: (1) First = FALSE Read dot/crs lon/lat
!              (2) Second = TRUE Read file
!

        if (lLALO) then
          if (ndim == 3) then
              sample = data( end_index(1),end_index(2),end_index(3),1 )
          else if (ndim == 2) then
              sample = data( end_index(1),end_index(2),1,1)
          else if (ndim == 1) then
              sample = data( end_index(1),1,1,1)
          end if

          if (lmore) write(*,'(A8,1x,I1,4(1x,I3),1x,A,1x,A," : ", F20.8,1x,A)')&
               name, ndim, end_index(1), end_index(2), end_index(3), end_index(4),&
               staggering, ordering, sample, trim(units)
!
! Call write_grib, if it is necessary
!
           innn=0
           innn_end=0
           ltogribyes= .false.
           do inn=1, inn_end
             if (name.eq.nametogrib(inn)) then

                if (((name.eq.'U        ').or.(name.eq.'V        ').or.&
                    (name.eq.'U10      ').or.(name.eq.'V10      ')).and.&
                    (ltogrib(inn))) then
                  if (name.eq.'U        ') nameT='V        '
                  if (name.eq.'V        ') nameT='U        '
                  if (name.eq.'U10      ') nameT='V10      '
                  if (name.eq.'V10      ') nameT='U10      '
                  do innT=1, inn_end
                    if ((nameT.eq.nametogrib(innT)).and.&
                        (sksec(inn)%seven.eq.sksec(innT)%seven).and.&
                        (sksec(inn)%eight.eq.sksec(innT)%eight).and.&
                        (sksec(inn)%nine.eq.sksec(innT)%nine)) then
                       innn=innn+1
                       in(innn)=inn
                       ltogribyes=((ltogribyes).or.(ltogrib(inn)))
                    endif
                  enddo
                elseif ((name.eq.'RAIN NON ').or.(name.eq.'RAIN CON ')) then
                  innn=innn+1
                  in(innn)=inn
                  ltogribyes=((ltogribyes).or.(ltogrib(inn)))
                  if (name.eq.'RAIN NON ') lRainNon = ltogrib(inn) 
                  if (name.eq.'RAIN CON ') lRainCon = ltogrib(inn)
                  nameRT='RAIN TOT '
                  do innT=1, inn_end
                    if (nameRT.eq.nametogrib(innT)) then
                       rksecRT=sksec(innT)
                       lRainTot = ltogrib(innT)
                       ltogribyes=((ltogribyes).or.(ltogrib(innT)))
                    endif
                  enddo
                else
                  innn=innn+1
                  in(innn)=inn
                  ltogribyes=((ltogribyes).or.(ltogrib(inn)))
                endif

             endif
           enddo
           innn_end=innn
!           print*,"Hola pepe ", innn_end, " Hola pedro ",ltogribyes

           if (ltogribyes) then
             do innn=1, innn_end
                rksec(innn)=sksec(in(innn))
                if (lmore) write(*,'(/,8x, "MM5toGRIB: Variable ", A, 1x, "selected.")') nametogrib(in(innn))
             enddo
             if (lmore) write(*,'(10x,A8,1x,I1,4(1x,I3),1x,A,1x,A," : ", F20.8,1x,A,/)')&
                        name, ndim, end_index(1), end_index(2), end_index(3), end_index(4),&
                        staggering, ordering, sample, trim(units)
!             stop
           if ((lTRain.eq.0).and.(lRainTot)) then
              if (allocated(dataTRain)) then
                deallocate(dataTRain)
              endif
              allocate(dataTRain(end_index(1), end_index(2), 1, 1))
              dataTRain=0
           endif
           if (llonlat) then
              if (lweight) then
                call grid_lonlat_input(end_index_n(1),end_index_n(2),&
                                       start_lon,start_lat,interval_lon,interval_lat,inter,lmore)
                allocate(wiscrs(end_index_n(1),end_index_n(2),4))
                allocate(wisdot(end_index_n(1),end_index_n(2),4))
                allocate(wagdot(end_index_n(1),end_index_n(2),1))
                allocate(index(end_index_n(1),end_index_n(2),4,2))
                allocate(indexUV(end_index_n(1),end_index_n(2),4,2))
                allocate(datan(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4)))
                allocate(datanUV(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4), 2))
              endif ! lweight
               end_index_n(3)=end_index(3)
               end_index_n(4)=end_index(4)
               deallocate(datan)
               allocate(datan(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4)))
               if (lUV.eq.0) then
                deallocate(datanUV)
                allocate(datanUV(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4), 2))
               endif
              if ((name.eq.'U        ').or.(name.eq.'V        ').or.&
                  (name.eq.'U10      ').or.(name.eq.'V10      ')) then
                lcross= .false.
                lUV=lUV+1
                call interpol_lonlat(bhi(7,1),bhr,loncrs,latcrs,londot,latdot,data,datan,datanUV,&
                                     name,index,indexUV,wiscrs,wisdot,wagdot,lmore)
                if ((name.eq.'U        ').or.(name.eq.'U10      ')) rksecU=rksec
                if ((name.eq.'V        ').or.(name.eq.'V10      ')) rksecV=rksec
                if (lUV.eq.2) then
                   if ((name(2:3)).eq.'10') then
                   name='U10      '
                   call write_grib(10,datanUV(1:end_index_n(1),1:end_index_n(2),1:end_index_n(3),1,2),&
                                   end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                   rksecU,PRESSURE)
                   name='V10      '
                   call write_grib(10,datanUV(1:end_index_n(1),1:end_index_n(2),1:end_index_n(3),1,1),&
                                   end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                   rksecV,PRESSURE)
!                   print*, 'Hola U10;V10'
                   else
                   name='U        '
                   call write_grib(10,datanUV(1:end_index_n(1),1:end_index_n(2),1:end_index_n(3),1,2),&
                                   end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                   rksecU,PRESSURE)
                   name='V        '
                   call write_grib(10,datanUV(1:end_index_n(1),1:end_index_n(2),1:end_index_n(3),1,1),&
                                   end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                   rksecV,PRESSURE)
!                   print*, 'Hola U;V'
                   endif
                   lUV=0
                endif
              elseif ((name.eq.'RAIN NON ').or.(name.eq.'RAIN CON ')) then
                lcross= .true.
                lUV=0
                if ((lRainNon).or.(lRainCon)) then
                call interpol_lonlat(bhi(7,1),bhr,loncrs,latcrs,londot,latdot,data,datan,datanUV,&
                                     name,index,indexUV,wiscrs,wisdot,wagdot,lmore)
                call write_grib(10,datan,end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                rksec,PRESSURE)
                lRainNon = .false.
                lRainCon = .false.
                endif
                if (lRainTot) then
                  if (name.eq.'RAIN NON ') then
                    dataTRain = dataTRain + data
                    lTRain=lTRain+1
                  endif
                  if (name.eq.'RAIN CON ') then
                    dataTRain = dataTRain + data
                    lTRain=lTRain+3
                  endif
                  if (lTRain.eq.4) then
                    call interpol_lonlat(bhi(7,1),bhr,loncrs,latcrs,londot,latdot,dataTRain,datan,datanUV,&
                                         name,index,indexUV,wiscrs,wisdot,wagdot,lmore)
                    name=nameRT  ! 'RAIN TOT '
                    call write_grib(10,datan,end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                    rksecRT,PRESSURE)
                    lTRain=0
                  elseif ((lTRain.ne.1).and.(lTRain.ne.3)) then
                    lTRain=0
                  endif
                endif
              else
                lcross= .true.
                lUV=0
                if ((name.ne.'PRESSURE ').and.(name.ne.'LATITDOT ').and.(name.ne.'LONGIDOT ').and.&
                     (name.ne.'LATITCRS ').and.(name.ne.'LONGICRS ')) then
                call interpol_lonlat(bhi(7,1),bhr,loncrs,latcrs,londot,latdot,data,datan,datanUV,&
                                     name,index,indexUV,wiscrs,wisdot,wagdot,lmore)
                call write_grib(10,datan,end_index_n,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                rksec,PRESSURE)
                else
                call write_grib(10,data,end_index,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                rksec,PRESSURE)
                endif
              endif
!             call interpol_lonlatNAG
!              stop
           else ! Not llonlat
              if ((name.eq.'RAIN NON ').or.(name.eq.'RAIN CON ')) then
                if ((lRainNon).or.(lRainCon)) then
                call write_grib(bhi(7,1),data,end_index,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                rksec,PRESSURE)
                lRainNon = .false.
                lRainCon = .false.
                endif
                if (lRainTot) then
                  if (name.eq.'RAIN NON ') then
                    dataTRain = dataTRain + data
                    lTRain=lTRain+1
                  endif
                  if (name.eq.'RAIN CON ') then
                    dataTRain = dataTRain + data
                    lTRain=lTRain+3
                  endif
                  if (lTRain.eq.4) then
                    name=nameRT  ! 'RAIN TOT '
!                    print*,"Hola name",name
                    call write_grib(bhi(7,1),dataTRain,end_index,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                    rksecRT,PRESSURE)
                    lTRain=0
                  elseif ((lTRain.ne.1).and.(lTRain.ne.3)) then
                    lTRain=0
                  endif
                endif
              else
                call write_grib(bhi(7,1),data,end_index,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                                rksec,PRESSURE)
              endif
           endif
           else
              if (innn.eq.0) then
                 write(*,'(/,8x, "MM5toGRIB:  Variable ", A, 1x, "not found in namelist.output.",/)') name
!                 stop
              else
                 if (lmore) then
                   do innn=1, innn_end
                      write(*,'(/,8x, "MM5toGRIB: Variable ", A, 1x, "not selected."/)') nametogrib(in(innn))
                   enddo
                 endif 
              endif
           endif
        else
          select case (name)
            case ('LATITCRS ')
              LATITCRS=nint(1000.*data(1,1,1,1))
              ENDLATITCRS=nint(1000.*data(end_index(1),end_index(2),1,1))
              if (lweight) then
                allocate(latcrs(end_index(1), end_index(2)))
                latcrs(1:end_index(1),1:end_index(2))=data(1:end_index(1), 1:end_index(2), 1, 1)
              endif 
            case ('LONGICRS ')
              LONGICRS=nint(1000.*data(1,1,1,1))
              ENDLONGICRS=nint(1000.*data(end_index(1),end_index(2),1,1))
              if (lweight) then
                allocate(loncrs(end_index(1), end_index(2)))
                loncrs(1:end_index(1),1:end_index(2))=data(1:end_index(1), 1:end_index(2), 1, 1)
              endif
            case ('LATITDOT ')
              LATITDOT=nint(1000.*data(1,1,1,1))
              ENDLATITDOT=nint(1000.*data(end_index(1),end_index(2),1,1))
              if (lweight) then
                allocate(latdot(end_index(1), end_index(2)))
                latdot(1:end_index(1),1:end_index(2))=data(1:end_index(1), 1:end_index(2), 1, 1)
              endif
            case ('LONGIDOT ')
              LONGIDOT=nint(1000.*data(1,1,1,1))
              ENDLONGIDOT=nint(1000.*data(end_index(1),end_index(2),1,1))
              if (lweight) then
                allocate(londot(end_index(1), end_index(2)))
                londot(1:end_index(1),1:end_index(2))=data(1:end_index(1), 1:end_index(2), 1, 1)
              endif
            case ('PRESSURE ')
!              write(*,'("bhi(12,8)=",I13)') bhi(12,8)
              do ipp=1, bhi(12,8)
                 PRESSURE(ipp)=(nint(data(ipp,1,1,1)))/100
              enddo
          end select
          if ((LATITCRS.NE.99999).AND.(LONGICRS.NE.999999).AND.&
             &(LATITDOT.NE.99999).AND.(LONGIDOT.NE.999999).AND.&
             &(PRESSURE(1).NE.-99999)) then
              close(iunit, err=98)
              open(iunit, file=flnm, err=99, form='unformatted', status='old', action='read')
              lLALO= .TRUE.
!              print*, 'hola end_index', end_index(1), end_index(2)
!              print*, 'hola dot', londot(80,94), latdot(80,94)
!              print*, 'hola crs', loncrs(80,94), latcrs(80,94)
!              stop
          endif
        endif

        deallocate(data)

     elseif (flag == 2) then
        newtime = .TRUE.
     else
        if (.NOT.(lLALO)) then
          write(*,'(/,8x,"Not found LATITDOT, LATITCRS, LONGIDOT,  LONGICRS and PRESSURE&
                     & from MM5 file.",/,8x,"MM5toGRIB: Impossible to write a grib.")')
        endif 
        stop
     endif
     read(iunit, iostat=ierr) flag
  enddo

  goto 100

98  write(*,'(/,8x, "Cannot close the ", A, 1x, "file",/)') trim(flnm)
    call abort()
    stop

99  write(*,'(/,8x, "Cannot open the ", A, 1x, "file",/)') trim(flnm)
    call help
    stop

100  continue
     if ((lweight).and.(ierrr == 0)) then
        if (allocated(loncrs)) deallocate(loncrs)
        if (allocated(latcrs)) deallocate(latcrs)
        if (allocated(londot)) deallocate(londot)
        if (allocated(latdot)) deallocate(latdot)
     endif

     if ((llonlat).and.(ierrr == 0)) then
!        print*, 'hola'
        deallocate(wiscrs)
        deallocate(wisdot)
        deallocate(wagdot)
        deallocate(index)
        deallocate(indexUV)
        deallocate(datan)
        deallocate(datanUV)
     endif

     if (lmore) write(*,'(/,"Hit the end of file of unit ", I3)') iunit
     if (lmore) write(*,'(/,8x, "Successful end of MM5toGRIB",/)')
     close(iunit, err=98)
     if (ierrr == -1) then
        write(*,'(/,8x, "Warning: perhaphs the ", A, 1x, "file is empty",/)') trim(flnm)
        stop
     endif

!     stop

end program MM5toGRIB

   subroutine printout_big_header(bhi, bhr, bhic, bhrc)

     implicit none
     integer, dimension(50,20) :: bhi
     real, dimension(20,20) :: bhr
     character(len=80), dimension(50,20) :: bhic
     character(len=80), dimension(20,20) :: bhrc
     integer :: i, j, v3j
   
     write(*,'(/)')
     v3j = bhi(1,1)
     if (bhi(1,1) == 11) v3j = v3j+5
     do j = 1, v3j
      if (j < 8 .or. j>10) then
        if (j == 1) write(*, '("TERRAIN Portion of big header:")')
        if (j == 2) write(*, '(/,"REGRID Portion of big header:")')
        if (j == 3) write(*, '(/,"RAWINS Portion of big header:")')
        if (j == 4) write(*, '(/,"SFC RAWINS Portion of big header:")')
        if (j == 5) write(*, '(/,"INTERP Portion of big header:")')
        if (j == 11) write(*, '(/,"MM5 Portion of big header:")')
        if (j == 6) write(*, '(/,"MM5 Substrate Temp File big header:")')
        if (j == 7) write(*, '(/,"MM5 Boundary File big header:")')
        if (j == 8) write(*, '(/,"Interpolated MM5 Portion of big header:")')
        write(*,'(/,"***Integers:"/)')
        do i = 1, size(bhi,1)
           if (bhi(i,j) /= -999) then
              write(*,'("BHI(",I3,",",I2,"):",I8," : ",A)')&
                   i, j, bhi(i,j),trim(bhic(i,j))
           endif
        enddo
   
        write(*,'(/,"***Floats:"/)')
        do i = 1, size(bhr,1)
           if (bhr(i,j) /= -999.) then
              write(*,'("BHR(",I3,",",I2,"):",F9.2," : ",A)')&
                   i, j, bhr(i,j),trim(bhrc(i,j))
           endif
        enddo
        write(*,'(/)')
      endif
     enddo
   end subroutine printout_big_header
 
   subroutine arguments(v2file, lmore, lsplit, llonlat)
     implicit none
     character(len=*) :: v2file
     character(len=120) :: harg
     logical :: lmore
     logical :: lsplit
     logical :: llonlat
   
#if defined CRAY
     integer(kind=4) :: ierr, i, numarg, ilen, ierror
     integer(kind=4), external :: ipxfargc
#else
     integer :: ierr, i, numarg
     integer, external :: iargc
#endif

#if defined CRAY
     numarg = ipxfargc()
#else 
     numarg = iargc()
#endif

     i = 1 
!     print*, 'hola', i , numarg  

     if ((numarg == 0).or.(numarg > 2)) then
       call help
     endif

     lmore = .false.
     lsplit= .false.
     llonlat= .false.

#if defined CRAY
     call pxfgetarg(i, harg, ilen, ierror)
#else
     call getarg(i, harg)
#endif

     if (harg == "-v") then
        i = i + 1
        if (numarg == 2) then
           lmore = .true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
        else
           call help
        endif
     elseif (harg == "-l") then
        i = i + 1
        if (numarg == 2) then
           llonlat= .true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
        else
           call help
        endif
     elseif ((harg == "-vl").or.(harg == "-lv")) then
        i = i + 1
        if (numarg == 2) then
           lmore = .true.
           llonlat=.true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
        else
           call help
        endif
     elseif ((harg == "-vo").or.(harg == "-ov")) then
        i = i + 1
        if (numarg == 2) then
           lmore = .true.
           lsplit= .true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
        else
           call help
        endif
     elseif ((harg == "-vol").or.(harg == "-ovl").or.(harg == "-vlo").or.(harg == "-olv")&
         .or.(harg == "-lvo").or.(harg == "-lov")) then
         i = i + 1
         if (numarg == 2) then
            lmore = .true.
            lsplit= .true.
            llonlat= .true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
         else
            call help
         endif
     elseif (harg == "-h") then
           call help
     elseif (harg == "-o") then
        i = i + 1
        if (numarg == 2) then
           lsplit= .true.
#if defined CRAY
           call pxfgetarg(i, harg, ilen, ierror)
#else
           call getarg(i, harg)
#endif
        else
           call help
        endif
      elseif ((harg == "-lo").or.(harg == "-ol")) then
         i = i + 1
         if (numarg == 2) then
            lsplit= .true.
            llonlat= .true.
#if defined CRAY
            call pxfgetarg(i, harg, ilen, ierror)
#else
            call getarg(i, harg)
#endif
         else
            call help
         endif
     endif

!#if defined CRAY
!        call pxfgetarg(i, harg, ilen, ierror)
!#else
!        call getarg(i, harg)
!#endif
!        print*, 'harg = ', trim(harg)
!
!        if (harg == "-v") then
!           i = i + 1
!           lmore = .true.
!        elseif (harg == "-h") then
!           call help
!        endif
!   
!     enddo
!
!#if defined CRAY
!     call pxfgetarg(i, harg, ilen, ierror)
!#else
!     call getarg(i, harg)
!#endif

     v2file = harg
   end subroutine arguments
   
   subroutine help
     implicit none
     character(len=120) :: cmd
     call getarg(0, cmd)
   
     write(*,'(/,"Usage: ", A, " [-v] MM5file ")') trim(cmd)
     write(*,'(8x, "-v       : Print info")')
     write(*,'(8x, "-o [-vo] : Split the gribfile for each period of time")')
     write(*,'(8x, "-l [-lvo]: Interpolate to regular longitude-latitude grid with")')
     write(*,'(19x, "the parameters of grid_lonlat.input file")')
     write(*,'(8x, "-h       : Print this help message.")')
     write(*,'(8x, "MM5file  : MMOUTP (isobaric) file name to read.",/)')
     write(*,'("Information code sources:")')
     write(*,'(13x, "readv3.f (MM5 Utils)")')
     write(*,'(13x, "module_horiz_interp.F (from REGRID MM5)")')
     write(*,'(13x, "map_utils_module.F (from LITTLE_R MM5)")')
     write(*,'(13x, "subroutine GRIBEX (EMOS library) from ECMWF for WMO FM-92 GRIB code",/)')
     write(*,'("2004 Instituto Nacional de Meteorologia (Spain)")')
     write(*,'("Servicio de Predecibilidad y Predicciones Extendidas")')
     write(*,'("Author: Alfons Callado Pallares")')
     write(*,'("E-mail: acallado@inm.es",/)')

     stop
   end subroutine help

   subroutine write_grib(igridtype,mm5data,end_index,lsplit,current_date,loneg,llonlat,lmore,NAMEFILE,innn_end,&
                         srksec,PRESSURE)
      implicit none
      common /KSEC/ KSEC0,KSEC1,KSEC2,KSEC3,KSEC4,PSEC2,PSEC3,ip,factor
!
!     Arrays are dimensioned
!     Integer    K.
!     Real       P.
!     Logical    O.
!     Character  H.
!
      INTEGER JPACK, JPBYTES, Pbytes
!      PARAMETER (JPACK=280000)
      PARAMETER (JPACK=100000)
      PARAMETER (JPBYTES=4)
#if defined REAL8
      parameter (Pbytes=8)
#else
      parameter (Pbytes=4)
#endif
!
!     Character argument
!
      character*1 :: HOPER = 'C'
!
!     Integer arguments for GRIBEX call
!
      INTEGER KSEC0(2)
      INTEGER KSEC1(1024)
      INTEGER KSEC2(1024)
      INTEGER KSEC3(2)
      INTEGER KSEC4(512)
      integer :: ndim
      integer, dimension(4) :: end_index
      integer ip
      integer j,i
!
!     Floating point arguments for GRIBEX call
!
      REAL(Pbytes) :: PSEC2(512)
      REAL(Pbytes) :: PSEC3(2)
      REAL(Pbytes), allocatable, dimension(:) :: PSEC4
      real, dimension(end_index(1),end_index(2),end_index(3),end_index(4)) :: mm5data
      real factor
!
!     Array for the GRIB product
!
      INTEGER KGRIB
      DIMENSION KGRIB(JPACK)
!
      INTEGER FILE, JCOUNT, NUMERR, KLENG, LENOUT
      INTEGER KLENP, KRET, KWORD, LOOP
!
      integer :: igridtype
      logical :: lsplit
      logical :: loneg
      logical :: llonlat
      logical :: lmore
      character*40 :: NAMEFILE
      character (len=24) :: current_date

      type sksecone
              integer :: six
              integer :: seven
              integer :: eight
              integer :: nine
              real    :: factor
      end type sksecone

      integer :: innn,innn_end
      type(sksecone),dimension(innn_end) :: srksec
      type(sksecone) :: rksec

      integer :: NPRESSURE
      parameter (NPRESSURE=50)
      integer, dimension(NPRESSURE) :: PRESSURE

!
!     Open the grib file
!
      if (lsplit) then
        NAMEFILE='grib2file'
        read(current_date(1:16),'(A16)') NAMEFILE(10:34)
      elseif (loneg) then
        NAMEFILE='gribfile'
        read(current_date(1:16),'(A16)') NAMEFILE(9:33)
        loneg= .false.
      endif
      if (lmore) write(*,'(/,8x,"Subroutine write_grib: writing to the grib file:",2x, A,/)') NAMEFILE
      CALL PBOPEN( FILE, NAMEFILE, 'a', KRET)
      IF( KRET.NE.0 ) THEN
        WRITE (*, *) 'Subroutine write_grib: Return code from PBOPEN = ',KRET
        CALL PBCLOSE(FILE, KRET)
        STOP 'Fault in PBOPEN'
      ENDIF
!
!     For each selection of same variable and levels of pressure, if you choose the variable
!
      do innn=1, innn_end
      rksec=srksec(innn)
      do ip=1, end_index(3)

      if ((rksec%seven.eq.0).or.((ip.eq.1).and.(rksec%seven.gt.100)).or.((ip.ne.1).and.(rksec%seven.eq.100))) then

      if ((rksec%eight.eq.0).or.(rksec%eight.eq.PRESSURE(ip)).or.((rksec%eight.gt.0).and.(rksec%eight.le.20))) then 

!
!     Make the header grib
!
      call header_grib(igridtype,lmore,rksec)
      if (KSEC1(6).NE.0) then
!
!     Data values from data(J,I,K,t) [MM5 format] to PSEC4(*) [GRIB format]
!
      allocate(PSEC4(KSEC2(2)*KSEC2(3)))
      if (llonlat) then
      do j=1, KSEC2(3)
        do i=1, KSEC2(2)
          PSEC4((j-1)*KSEC2(2) + i) = (mm5data(i,j,ip,1))*factor
        enddo
      enddo
      if (lmore) then
      write(*,'("KSEC1(6)=",I13)') KSEC1(6)
      write(*,'("KSEC1(7)=",I13)') KSEC1(7)
      write(*,'("KSEC1(8)=",I13)') KSEC1(8)
      write(*,'("KSEC1(9)=",I13)') KSEC1(9)
      write(*,'("factor   ",F13.8)') factor
      write(*,'("PSEC4(",I13,")=",F20.8)') KSEC2(2)*KSEC2(3), PSEC4(KSEC4(1))
      write(*,'("mm5data(",I3,",",I3,",1,1)=",F20.8)') KSEC2(2), KSEC2(3),&
                                                      mm5data(KSEC2(2),KSEC2(3),ip,1)
      write(*,'("PSEC4(",I13,")=",F20.8)') 1,  PSEC4(1)
      write(*,'("mm5data(",I3,",",I3,",1,1)=",F20.8)') 1,1,mm5data(1,1,ip,1)
      endif
      else
      do j=1, KSEC2(3)
        do i=1, KSEC2(2)
          PSEC4((j-1)*KSEC2(2) + i) = (mm5data(j,i,ip,1))*factor
        enddo
      enddo
      if (lmore) then
      write(*,'("KSEC1(6)=",I13)') KSEC1(6)
      write(*,'("KSEC1(7)=",I13)') KSEC1(7)
      write(*,'("KSEC1(8)=",I13)') KSEC1(8)
      write(*,'("Factor  =",F13.8)') factor
      write(*,'("PSEC4(",I13,")=",F20.8)') KSEC2(2)*KSEC2(3), PSEC4(KSEC4(1))
      write(*,'("mm5data(",I3,",",I3,",1,1)=",F20.8)') KSEC2(3), KSEC2(2),&
                                                      mm5data(KSEC2(3),KSEC2(2),ip,1)
      write(*,'("PSEC4(",I13,")=",F20.8)') 1,  PSEC4(1)
      write(*,'("mm5data(",I3,",",I3,",1,1)=",F20.8)') 1,1,mm5data(1,1,ip,1)
      endif
      endif
!
!       Code (pack) GRIB.
!
!      KLENG = (LENOUT+JPBYTES-1)/JPBYTES
       KLENG = JPACK
!      KLENP = JPACK*5
       KLENP = KSEC4(1)
!       Return to the calling routine even if an error is encountered
       KRET = 1
!      KRET = 0
      if (lmore) print*, 'Call GRIBEX'
      CALL GRIBEX(KSEC0,KSEC1,KSEC2,PSEC2,KSEC3,PSEC3,KSEC4,&
                  PSEC4,KLENP,KGRIB,KLENG,KWORD,HOPER,KRET)
      deallocate(PSEC4)
!
!      Check return code.
!
      IF (KRET .LT. 0 ) THEN
         WRITE (*, *) 'Warning subroutine write_grib: Return code from GRIBEX = ',KRET
      ELSEIF (KRET .GT. 0 ) THEN
          WRITE (*, *) 'Error subroutine write_grib: Return code from GRIBEX = ',KRET 
         CALL PBCLOSE(FILE, KRET)
         STOP 'Fault in GRIBEX'
      ELSEIF (lmore) THEN
         WRITE (*,*) 'Subroutine write_grib: GRIBEX return code = ',KRET 
      ENDIF
      IF (lmore) THEN
      WRITE (*,*) 'Subroutine write_grib:'
      WRITE (*,*) 'KLENP -> The number of elements in array PSEC4 = ',KLENP
      WRITE (*,*) 'KLENG -> The number of elements in array kGRIB = ',KLENG
      WRITE (*,*) 'KWORD -> Tnfe  of KGRIB occupied by coded data = ',KWORD
      ENDIF
!
!     Write the grib file
!
!      CALL PBWRITE(FILE,KGRIB,KWORD*8,KRET)
!      CALL PBWRITE(FILE,KGRIB,KWORD*16,KRET)
!      CALL PBWRITE(FILE,KGRIB,(KWORD*JPBYTES),KRET)
      CALL PBWRITE(FILE,KGRIB,(KWORD*JPBYTES)-JPBYTES+1,KRET)
!      CALL PBWRITE(FILE,KGRIB,(KWORD*JPBYTES)-2*JPBYTES,KRET)
!      CALL PBWRITE(FILE,KGRIB,KWORD*3,KRET)
      IF ((KRET .GE. 0 ).AND.(lmore)) THEN
         WRITE (*,*) 'Subroutine write_grib:'
         WRITE (*,*) '  PBWRITE return number of bytes written = ',KRET
      ELSEIF (KRET .EQ. -1) THEN
         WRITE (*,*) 'Subroutine write_grib:'
         WRITE (*,*) '  PBWRITE error writing to the file',KRET
         CALL PBCLOSE(FILE, KRET)
         STOP 'Fault in PBWRITE'
      ENDIF
      endif

      endif
      endif
      enddo
      enddo
!
!     Close the grib file
!
      CALL PBCLOSE(FILE, KRET)
!      stop
   end subroutine write_grib

   subroutine header_grib(igridtype,lmore,rksec)
      implicit none
      common /MM5SEC/ ordering, name, staggering, current_date, bhi, time,&
                      PRESSURE, LATITDOT, LATITCRS, LONGIDOT, LONGICRS,&
                      ENDLATITDOT, ENDLATITCRS, ENDLONGIDOT, ENDLONGICRS, bhr,&
                      nfield
      common /interpol/ lweight,lcross,lUV,end_index,end_index_n,&
                        inter,start_lon,start_lat,interval_lon,interval_lat
      common /KSEC/ KSEC0,KSEC1,KSEC2,KSEC3,KSEC4,PSEC2,PSEC3,ip,factor
      integer :: igridtype
      character (len= 4) ordering
      character (len= 9) name
      character (len= 4) staggering
      character (len=24) current_date
      integer ip
      integer :: NPRESSURE
      parameter (NPRESSURE=50)
      integer, dimension(NPRESSURE) :: PRESSURE
      integer, dimension(50,20) :: bhi
      integer :: LATITDOT,LATITCRS,LONGIDOT,LONGICRS
      integer :: ENDLATITDOT,ENDLATITCRS,ENDLONGIDOT,ENDLONGICRS
      real, dimension(20,20) :: bhr
      real :: time
      logical :: lweight, lcross, lUV, lmore
      integer, dimension(4) :: end_index,end_index_n
      real :: start_lon,start_lat,interval_lon,interval_lat
      integer :: inter
      integer il
      INTEGER KSEC0(2)
      INTEGER KSEC1(1024)
      INTEGER KSEC2(1024)
      INTEGER KSEC3(2)
      INTEGER KSEC4(512)
      integer :: Pbytes
#if defined REAL8
      parameter (Pbytes=8)
#else
      parameter (Pbytes=4)
#endif
      REAL(Pbytes) :: PSEC2(512)
      REAL(Pbytes) :: PSEC3(2)
      type sksecone
              integer :: six
              integer :: seven
              integer :: eight
              integer :: nine
              real    :: factor
      end type sksecone
      type(sksecone) :: rksec
      integer :: nfield
      real factor 
!         Inicialization
       data KSEC0,KSEC3/4*0/,KSEC1,KSEC2/2048*0/,KSEC4/512*0/
      if ((ordering.EQ.'YXP ').OR.(ordering.EQ.'YX  ').OR.(ordering.EQ.'P   ')) then
!         SECTION 0: THE INDICATOR SECTION
!         Number of octets in the GRIB message 
!      KSEC0(1) = 0
!         GRIB edition number
      KSEC0(2)=1
!         SECTION 1: THE PRODUCT DEFINITION SECTION (PDS)
!         Version number of code table 2 (=1 Pressure [Pa])
      KSEC1(1)=1
!         Identification of centre (=96 Athens)
      KSEC1(2)=96
!         Generating process identification number, allocated by the originating centre
      KSEC1(3)=1
!         Grid definition. NNN catalogue number of grid used by the originating centre
!         255 for a non-catalogued grid, in which case the grid description follows in KSEC2. 
      KSEC1(4)=255
!         Flag showing whether sections 2 and 3 are present (see WMO code table 1):
!         0,   Sections 2 and 3 are omitted.
!         64,  Section 2 is omitted, Section 3 is included.
!         128, Section 2 is included, Section 3 is omitted (bitmap section).
!         192, Sections 2 and 3 are included.
      KSEC1(5)=128
!      ip=2
      KSEC1(6)=rksec%six
      KSEC1(7)=rksec%seven
      KSEC1(8)=rksec%eight
      KSEC1(9)=rksec%nine
      factor=rksec%factor
      if (rksec%seven.eq.0) then
      select case (ip)
        case (1)
          KSEC1(7)=105
        case (2:NPRESSURE)
          KSEC1(7)=100
          KSEC1(8)=PRESSURE(ip)
        case default
          stop
      end select
      endif
!         Parameter indicator / Type of level indicator / Height, pressure, etc of level
      if (rksec%factor.eq.0) factor=1.0
      if (rksec%six.ne.0) KSEC1(6)=rksec%six
!        Year of century (YY) / Month (MM) / Day (DD) / Hour (HH) / Minute (MM)
!        Time unit indicator (KSEC1(15)=0 minute; =1 hour; etc)
!        Century of reference time of data (eg 20 for 1997)
!      read(current_date(3:4),'(I2)') KSEC1(10)
!      read(current_date(6:7),'(I2)') KSEC1(11)
!      read(current_date(9:10),'(I2)') KSEC1(12)
!      read(current_date(12:13),'(I2)') KSEC1(13)
!      read(current_date(15:16),'(I2)') KSEC1(14)
!     Analysis time (of forecast)
      KSEC1(10)=BHI(5,11)-((int(BHI(5,11)/100))*100)
      KSEC1(11)=BHI(6,11)
      KSEC1(12)=BHI(7,11)
      KSEC1(13)=BHI(8,11)
      KSEC1(14)=BHI(9,11)
!     Unit of time = Hour
      KSEC1(15)=1
!     Forecast time
      KSEC1(16)=anint(time/60.)
      KSEC1(17)=0
      KSEC1(18)=0
!      read(current_date(1:2),'(I2)') KSEC1(21)
      KSEC1(21)= int(BHI(5,11)/100)
      KSEC1(21)= 1 + KSEC1(21)
!      print*,"Hola select name",name
      select case (name)
        case ('U        ')  
          if (rksec%six.eq.0) KSEC1(6)=34
        case ('V        ')
          if (rksec%six.eq.0) KSEC1(6)=33
        case ('T        ')
          if (rksec%six.eq.0) KSEC1(6)=11
        case ('Q        ')
          if (rksec%six.eq.0) KSEC1(6)=53
!        Cloud mixing ratio kg/kg , WMO 76 --> Cloud water Kg/m**2
        case ('CLW      ')
          if (rksec%six.eq.0) KSEC1(6)=128
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Cloud mixing ratio kg/kg , not WMO 76 --> Cloud water Kg/m**2.",/)') name,KSEC1(6)
!        Rain water mixing ratio Kg/Kg, WMO 63 --> Total precipitation Kg/m**2
        case ('RNW      ')
          if (rksec%six.eq.0) KSEC1(6)=129
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Rain water mixing ratio Kg/Kg, not WMO 63 --> Total precipitation Kg/m**2.",/)') name,KSEC1(6)
!        Atmospheric radiative tendency K/day
        case ('RAD TEND ')
          if (rksec%six.eq.0) KSEC1(6)=130
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Atmospheric radiative tendency K/day.",/)') name,KSEC1(6)
!        Vertical wind component m/s, WMO 40 --> Vertical velocity m/s
        case ('W        ')
          if (rksec%six.eq.0) KSEC1(6)=40
!        Pressure perturbation Pa, WMO 26 ? --> Pressure anomaly Pa 
        case ('PP       ')
          if (rksec%six.eq.0) KSEC1(6)=131
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Pressure perturbation Pa, WMO 26 ? --> Pressure anomaly Pa.",/)') name,KSEC1(6)
!        Altitude geopotencial --> Geopotencial ?
        case ('H        ')
          if (rksec%six.eq.0) then
             KSEC1(6)=6
             factor=9.80665
          endif
!        Relative humidity
        case ('RH       ')
          if (rksec%six.eq.0) KSEC1(6)=52
!          factor=0.01
!        (REFERENCE) SURFACE PRESSURE MINUS PTOP
        case ('PSTARCRS ')
          if (rksec%six.eq.0) KSEC1(6)=132
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," (REFERENCE) SURFACE PRESSURE MINUS PTOP.",/)') name,KSEC1(6)
!        Ground temperature K
        case ('GROUND T ')
          if (rksec%six.eq.0) KSEC1(6)=11
          if (rksec%eight.eq.0) KSEC1(8)=0
!        Accumulated Convective Precipition cm, WMO 63 --> Convective Precipitation Kg/m**2
        case ('RAIN CON ')
!          KSEC1(6)=133
          if (rksec%six.eq.0) then
             KSEC1(6)=63
             factor=10.
          endif
          KSEC1(16)=0
          KSEC1(17)=anint(time/60.)
          KSEC1(18)=4
!          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
!                   &"is 133 Accumulated Convective Precipition cm, WMO 63 --> Convec. Precip. Kg/m**2.",/)') name
!        Accumulated nonConvective Precipition cm
        case ('RAIN NON ')
!          KSEC1(6)=134
          if (rksec%six.eq.0) then
             KSEC1(6)=62
             factor=10.
          endif
          KSEC1(16)=0
          KSEC1(17)=anint(time/60.)
          KSEC1(18)=4
!          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
!                   &"is 134 Accumulated nonConvective Precipition cm.",/)') name
!       Total Accumulated Precipition cm (nonConvective + Convective)
        case ('RAIN TOT ')
          if (rksec%six.eq.0) then
             KSEC1(6)=61
             factor=10.
          endif
          KSEC1(16)=0
          KSEC1(17)=anint(time/60.)
          KSEC1(18)=4
!        [Cloud ice Kg/m**2, WMO 58 -->] Ice cloud mixing ratio (kg/kg) (if IMPHYS5)
        case ('ICE      ')
          if (rksec%six.eq.0) KSEC1(6)=160
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Ice cloud mixing ratio (kg/kg).",/)') name,KSEC1(6)
!        [Convective snow / Large scale snow Kg/m**2, WMO 78 / 79 -->] Snow mixing ratio (kg/kg) (if IMPHYS5)
        case ('SNOW     ')
          if (rksec%six.eq.0) KSEC1(6)=161
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Snow mixing ratio (kg/kg).",/)') name,KSEC1(6)
!        Graupel (kg/kg) (if IMPHYS6)
        case ('GRAUPEL  ')
          if (rksec%six.eq.0) KSEC1(6)=162
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Graupel (kg/kg).",/)') name,KSEC1(6)
!        Number concentration of ice (if IMPHYS=7)
        case ('NCI      ')
          if (rksec%six.eq.0) KSEC1(6)=163
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Number concentration of ice.",/)') name,KSEC1(6)
!        Terrain Elevation m
        case ('TERRAIN  ')
          if (rksec%six.eq.0) KSEC1(6)=135
          if (rksec%eight.eq.0) KSEC1(8)=0
!          factor=9.80665
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Terrain Elevation m.",/)') name,KSEC1(6)
!        MAP SCALE FACTOR (DIMENSIONLESS)
        case ('MAPFACCR ')
          if (rksec%six.eq.0) KSEC1(6)=136
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," MAP SCALE FACTOR (cross points) (DIMENSIONLESS).",/)') name,KSEC1(6)
!        MAP SCALE FACTOR (DIMENSIONLESS)
        case ('MAPFACDT ')
          if (rksec%six.eq.0) KSEC1(6)=137
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," MAP SCALE FACTOR (dot points) (DIMENSIONLESS).",/)') name,KSEC1(6)
!        CORIOLIS PARAMETER 1/s
        case ('CORIOLIS ')
          if (rksec%six.eq.0) KSEC1(6)=138
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," CORIOLIS PARAMETER 1/s.",/)') name,KSEC1(6)
!        INFINITE RESERVOIR SLAB TEMPERATURE k
        case ('RES TEMP ')
          if (rksec%six.eq.0) KSEC1(6)=139
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," INFINITE RESERVOIR SLAB TEMPERATURE k.",/)') name,KSEC1(6)
!        Latitude cross points (south negative) (DEGREES)
        case ('LATITCRS ')
          if (rksec%six.eq.0) KSEC1(6)=140
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Latitude cross points (south negative) (DEGREES).",/)') name,KSEC1(6)
!        Longitude cross points (west negative) (DEGREES)
        case ('LONGICRS ')
          if (rksec%six.eq.0) KSEC1(6)=141
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," Longitude cross points (west negative) (DEGREES).",/)') name,KSEC1(6)
!        LANDUSE CATEGORY (category)
        case ('LAND USE ')
          if (rksec%six.eq.0) KSEC1(6)=142
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," LANDUSE CATEGORY (category).",/)') name,KSEC1(6)
!        Sea surface temperature
        case ('TSEASFC  ')
          if (rksec%six.eq.0) KSEC1(6)=11
          if (rksec%seven.eq.0) KSEC1(7)=102
          if (rksec%eight.eq.0) KSEC1(8)=0
!        PBL height m, WMO 67 --> Mixed layer depth m
        case ('PBL HGT  ')
          if (rksec%six.eq.0) KSEC1(6)=67
          if (rksec%eight.eq.0) KSEC1(8)=0
!        PBL REGIME (DIMENSIONLESS)
        case ('REGIME   ')
          if (rksec%six.eq.0) KSEC1(6)=143
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," PBL REGIME (DIMENSIONLESS).",/)') name,KSEC1(6)
!        Sensible heat Flux W/m**2
        case ('SHFLUX   ')
          if (rksec%six.eq.0) KSEC1(6)=122
          if (rksec%eight.eq.0) KSEC1(8)=0
!        Latent heat flux W/m^2
        case ('LHFLUX   ')
          if (rksec%six.eq.0) KSEC1(6)=121
          if (rksec%eight.eq.0) KSEC1(8)=0
!        FRICTIONAL VELOCITY m/s
        case ('UST      ')
          if (rksec%six.eq.0) KSEC1(6)=144
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," FRICTIONAL VELOCITY m/s.",/)') name,KSEC1(6)
!        SURFACE DOWNWARD SHORTWAVE RADIATION W/m^2
        case ('SWDOWN   ')
          if (rksec%six.eq.0) KSEC1(6)=145
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," SURFACE DOWNWARD SHORTWAVE RADIATION W/m^2.",/)') name,KSEC1(6)
!        SURFACE DOWNWARD LONGWAVE RADIATION W/m^2
        case ('LWDOWN   ')
          if (rksec%six.eq.0) KSEC1(6)=146
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," SURFACE DOWNWARD LONGWAVE RADIATION W/m^2.",/)') name,KSEC1(6)
!        TOP OUTGOING SHORTWAVE RADIATION W/m^2
        case ('SWOUT    ')
          if (rksec%six.eq.0) KSEC1(6)=147
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," TOP OUTGOING SHORTWAVE RADIATION W/m^2.",/)') name,KSEC1(6)
!        TOP OUTGOING LONGWAVE RADIATION W/m^2
        case ('LWOUT    ')
          if (rksec%six.eq.0) KSEC1(6)=148
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4," TOP OUTGOING LONGWAVE RADIATION W/m^2.",/)') name,KSEC1(6)
!        SOIL TEMPERATURE IN LAYER  N
        case ('SOIL T 1 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=1
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 1.",/)') name,KSEC1(6),KSEC1(8)
        case ('SOIL T 2 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=2
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 2.",/)') name,KSEC1(6),KSEC1(8)
        case ('SOIL T 3 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=3
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 3.",/)') name,KSEC1(6),KSEC1(8)
        case ('SOIL T 4 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=4
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 4.",/)') name,KSEC1(6),KSEC1(8)
        case ('SOIL T 5 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=5
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 5.",/)') name,KSEC1(6),KSEC1(8)
        case ('SOIL T 6 ')
          if (rksec%six.eq.0) KSEC1(6)=149
          if (rksec%eight.eq.0) KSEC1(8)=6
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                   &"is ",I4,"/",I4," SOIL TEMPERATURE IN LAYER 6.",/)') name,KSEC1(6),KSEC1(8)
!        2-meter temperature, WMO 11, 105, 2
        case ('T2       ')
          if (rksec%six.eq.0) KSEC1(6)=11
          if (rksec%eight.eq.0) KSEC1(8)=2
!        2-meter mixing ratio
        case ('Q2       ')
          if (rksec%six.eq.0) KSEC1(6)=53
          if (rksec%eight.eq.0) KSEC1(8)=2
!        10-meter U Component
        case ('U10      ')
          if (rksec%six.eq.0) KSEC1(6)=34
          if (rksec%eight.eq.0) KSEC1(8)=10
!        10-meter V Component
        case ('V10      ')
          if (rksec%six.eq.0) KSEC1(6)=33
          if (rksec%eight.eq.0) KSEC1(8)=10
!        Pressure surface
        case ('PSFC     ')
         if (rksec%six.eq.0) KSEC1(6)=1
         if (rksec%eight.eq.0) KSEC1(8)=0
!        Pressure sea level in cross points
        case ('PSEALVLC ')
         if (rksec%six.eq.0) KSEC1(6)=1
         if (rksec%seven.eq.0) KSEC1(7)=103
         if (rksec%eight.eq.0) KSEC1(8)=0
!        Pressure sea level in dot points
        case ('PSEALVLD ')
         if (rksec%six.eq.0) KSEC1(6)=1
         if (rksec%seven.eq.0) KSEC1(7)=103
         if (rksec%eight.eq.0) KSEC1(8)=1
         if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4,"/",I4,"/",I4," Pressure sea level in dot points.",/)') name,KSEC1(6),KSEC1(7),KSEC1(8)
!        Latitude dot points (DEGREES)
        case ('LATITDOT ')
          if (rksec%six.eq.0) KSEC1(6)=150
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4," Latitude dot points (DEGREES).",/)') name,KSEC1(6)
!        Longitude dot points (DEGREES)
        case ('LONGIDOT ')
          if (rksec%six.eq.0) KSEC1(6)=151
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4," Longitude dot points (DEGREES).",/)') name,KSEC1(6)
!        Surface relative humidity
        case ('RH SFC   ')
          if (rksec%six.eq.0) KSEC1(6)=52
          if (rksec%seven.eq.0) KSEC1(7)=105
          if (rksec%eight.eq.0) KSEC1(8)=0
!        Pressure levels in grid
        case ('PRESSURE ')
          if (rksec%six.eq.0) KSEC1(6)=152
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4," Pressure levels in grid.",/)') name,KSEC1(6)
        case ('SUMFB    ')
          if (rksec%six.eq.0) KSEC1(6)=153
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4," Mass flux updraft (if ICUPA=8).",/)') name,KSEC1(6)
        case ('SPSRC    ')
          if (rksec%six.eq.0) KSEC1(6)=154
          if (lmore) write(*,'(/,8x, "MM5toGRIB / write_grib / header_grib:  warning, variable ",/,8x,A, 1x,&
                  &"is ",I4," Source layer updraft (if ICUPA=8).",/)') name,KSEC1(6)
        case default
          KSEC1(6)=rksec%six
          if (rksec%six.eq.0) then
            nfield=nfield+164
            KSEC1(6)=nfield
          endif
          if (rksec%seven.ne.0) then
            KSEC1(7)=rksec%seven
            KSEC1(8)=rksec%eight
          endif
!           KSEC1(6)=0 
!           return
      end select
!         SECTION 2: GRID DESCRIPTION SECTION (GDS)
!         Data representation type (see WMO code table 6):
!         5 for polar stereographic              --> case(2)
!         3 Lambert conformal                    --> case(1)
!         1 Mercator                             --> case(3)
!         0 Latitude/longitude grid (10 Rotated) --> case(10)
      select case (igridtype)
        case (1)
          KSEC2(1)=3
          KSEC2(2)=end_index(2)
          KSEC2(3)=end_index(1)
          if (staggering.EQ.'D   ') then
          KSEC2(4)=LATITDOT
          KSEC2(5)=LONGIDOT
          elseif (staggering.EQ.'C   ') then
          KSEC2(4)=LATITCRS
          KSEC2(5)=LONGICRS
          endif
!         Resolution flag:
!         0     Direction increments not given.
!         128   Direction increments given.
          KSEC2(7)=nint(1000.*bhr(3,1))
          KSEC2(9)=nint(bhr(9,1))
          KSEC2(10)=nint(bhr(9,1))
          KSEC2(11)=64
!         Number of vertical coordinate parameters ??
          KSEC2(12)=0
          KSEC2(13)=0
!         First latitude from the pole at which the secant cone cuts the sphere.
!         Second latitude at which the secant cone cuts the sphere.
          KSEC2(14)=nint(1000.*bhr(5,1))
          KSEC2(15)=nint(1000.*bhr(6,1))
          KSEC2(17)=0
          KSEC2(18)=0
          KSEC2(19)=0
          PSEC2(1)=0.
          PSEC2(2)=bhr(4,1)
        case (2)
          KSEC2(1)=5
!         Number of points along X-axis. bhi(6,1) or  bhi(10,1)
!         Number of points along Y-axis. bhi(5,1) or  bhi(9,1) 
          KSEC2(2)=end_index(2)
          KSEC2(3)=end_index(1)
!         Latitude of first grid point.
!         Longitude of first grid point.
          if (staggering.EQ.'D   ') then
          KSEC2(4)=LATITDOT
          KSEC2(5)=LONGIDOT
          elseif (staggering.EQ.'C   ') then
          KSEC2(4)=LATITCRS
          KSEC2(5)=LONGICRS
          endif
!         Orientation of the grid, ie the longitude of the meridian which is parallel
!         to the Y-axis along which latitude increases as the Y-coordinate increases.
          KSEC2(7)=nint(1000.*bhr(3,1))
!         X-direction grid length. bhr(1,1) or bhr(9,1) ?
!         bhr(9,1) includes both coarse domain and sub-domains
!         Y-direction grid length. (in meters)
          KSEC2(9)=nint(bhr(9,1))
          KSEC2(10)=nint(bhr(9,1))
!         Scanning mode flag (see WMO code table 8).
          KSEC2(11)=64
!         Projection centre flag:
!          0      North pole is on projection plane.
!          128    South pole is on projection plane.
          KSEC2(13)=0
!         0 Regular grid.
          KSEC2(17)=0
!         Earth flag:
!          0      Earth assumed spherical with radius 6367.47 km.
!          64     Earth assumed oblate spheroidal with size as determined by IAU in 1965:
!                    6378.160km,   6356.775km,   f=1/297.0
          KSEC2(18)=0
!         Components flag. Valid values are:
!          0      Resolved u and v components of vector quantities relative to easterly
!                 and northerly directions.
!          8      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. 
          KSEC2(19)=0
!         Angle of rotation
          PSEC2(1)=0.
!         Streching factor
          PSEC2(2)=bhr(4,1)
        case(3)
          KSEC2(1)=1
!         Ni - Number of points along a parallel:  end_index_n(1) E-W
!         Nj - Number of points along a meridian:  end_index_n(2) N-S
          KSEC2(2)=end_index(2)
          KSEC2(3)=end_index(1)
!         La1 - Latitude of first grid point.
!         Lo1 - Longitude of first grid point.
          if (staggering.EQ.'D   ') then
          KSEC2(4)=LATITDOT
          KSEC2(5)=LONGIDOT
          elseif (staggering.EQ.'C   ') then
          KSEC2(4)=LATITCRS
          KSEC2(5)=LONGICRS
          endif
!         Resolution flag:
!         0     Direction increments not given.
!         128   Direction increments given.
          KSEC2(6)=0
!         La1 - Latitude of last grid point.
!         Lo1 - Longitude of last grid point.
          if (staggering.EQ.'D   ') then
          KSEC2(7)=ENDLATITDOT
          KSEC2(8)=ENDLONGIDOT
          elseif (staggering.EQ.'C   ') then
          KSEC2(7)=ENDLATITCRS
          KSEC2(8)=ENDLONGICRS
          endif
!         Latitude at which the Mercator projection cylinder intersects the Earth(Latin).
          KSEC2(9)=nint(1000.*bhr(5,1))
!         Reserved. Set to 0.
          KSEC2(10)=0
!         Scanning mode flag (see WMO code table 8).
          KSEC2(11)=64
!         Number of vertical coordinate parameters ??
          KSEC2(12)=0
!         i direction grid length.
!         j direction grid length.
          KSEC2(13)=nint(bhr(9,1))
          KSEC2(14)=nint(bhr(9,1))
!         Reserved. Set to 0.
          KSEC2(15)=0
          KSEC2(16)=0
!         0 Regular grid.
          KSEC2(17)=0
!          0      Earth assumed spherical with radius 6367.47 km.
!          64     Earth assumed oblate spheroidal with size as determined by IAU in 1965:
!                    6378.160km,   6356.775km,   f=1/297.0
          KSEC2(18)=0
!         Components flag. Valid values are:
!          0      Resolved u and v components of vector quantities relative to easterly
!                 and northerly directions.
!          8      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.
          KSEC2(19)=0
!         Reserved. Set to 0.
          KSEC2(20)=0
          KSEC2(21)=0
          KSEC2(22)=0
        case(10)
          KSEC2(1)=0
!         Ni - Number of points along a parallel:  end_index_n(1) E-W  
!         Nj - Number of points along a meridian:  end_index_n(2) N-S
          KSEC2(2)=end_index_n(1)
          KSEC2(3)=end_index_n(2)
!         La1 - Latitude of first grid point.
!         Lo1 - Longitude of first grid point.
          KSEC2(4)=nint(1000.*start_lat)
          KSEC2(5)=nint(1000.*start_lon)
!         Resolution flag (0   --> Direction increments not given (quasi-regular grids)
!                         (128 --> Direction increments given (regular grids)
          KSEC2(6)=128
!         La2 - Latitude of last grid point.
!         Lo2 - Longitude of last grid point.
          KSEC2(7)=nint(1000.*(start_lat+(interval_lat*real((end_index_n(2))-1))))
          KSEC2(8)=nint(1000.*(start_lon+(interval_lon*real((end_index_n(1))-1))))
!         Di - i direction increment.
!         Dj - j direction increment.
          KSEC2(9)=nint(1000.*interval_lon)
          KSEC2(10)=nint(1000.*interval_lat)
!         Scanning mode flags (Code Table 8).
          KSEC2(11)=64
!         Number of vertical coordinate parameters. ??
          KSEC2(12)=0
!         Latitude of the southern pole of rotation. ??
          KSEC2(13)=0
!         Longitude of the southern pole of rotation. ??
          KSEC2(14)=0
!         Latitude of the pole of stretching. ??
          KSEC2(15)=0
!         Longitude of the pole of stretching. ??
          KSEC2(16)=0
!         0  Regular grid.
!         1  Quasi-regular (reduced) grid.
          KSEC2(17)=0
!         Earth flag:
!          0      Earth assumed spherical with radius 6367.47 km.
!          64     Earth assumed oblate spheroidal with size as determined by IAU in 1965:
!                    6378.160km,   6356.775km,   f=1/297.0
          KSEC2(18)=0
!         Components u,v flags:
!          0  Resolved u and v components of vector quantities relative to easterly
!             and northerly directions.
!          8  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.
          KSEC2(19)=0
        case default
          write(*,'(/,8x,"Subroutine write_grib / header_grib:",/, 20x,&
                        &"Error MM5toGRIB: Not Polar Stereograpic or Lambert Conformal or Mercator projection",/)')
          stop
      end select
!         SECTION 4: BINARY DATA SECTION (BDS)
!         Number of data values in array PSEC4 to be encoded.
      KSEC4(1)=KSEC2(2)*KSEC2(3)
!         Number of bits used for each packed value.
      KSEC4(2)=16
!         Type of packing:
!           0       Simple packing.
!           64      Complex or second-order packing.
      KSEC4(4)=0
!         Data representation:
!           0       Floating point data.
!           32      Integer data.
      KSEC4(5)=0
!         Additional flags indicator:
!           0       No additional flags.
!           16      Additional flags.
      KSEC4(6)=0
!         Number of values indicator:
!           0      Single datum at each grid point.
!           64     Matrix of values at each grid point.
      KSEC4(8)=0
!         Values width indicator:
!           0      Second order values have constant width.
!           16     Second order values have different widths.
      KSEC4(10)=0
!         Number of bits for second order values when these have constant width.
      KSEC4(11)=0
!      do il=1, 2
!         write(*,'("KSEC0(",I2,")=",I7)') il, KSEC0(il)
!      enddo
!      do il=1, 21
!        write(*,'("KSEC1(",I2,")=",I7)') il, KSEC1(il)
!      enddo
!      do il=1, 20
!         write(*,'("KSEC2(",I2,")=",I7)') il, KSEC2(il)
!      enddo
!      do il=1, 2
!         write(*,'("PSEC2(",I2,")=",F7.2)') il, PSEC2(il)
!      enddo
!      do il=1, 11
!          write(*,'("KSEC4(",I2,")=",I7)') il, KSEC4(il)
!      enddo
!      stop
      if (lmore) then
      write(*,'(/,8x,"MM5toGRIB:subroutine write_grib:subroutine header_grib: type of data and date",/)')
      write(*,'(10x,"KSEC0(",I2,")=",I7)') 1, KSEC0(1)
      write(*,'(10x,"KSEC0(",I2,")=",I7)') 2, KSEC0(2)
      write(*,'(10x,"KSEC2(",I2,")=",I7)') 6, KSEC2(6)
      write(*,'(10x,"KSEC2(",I2,")=",I7)') 7, KSEC2(7)
      write(*,'(10x,"KSEC2(",I2,")=",I7)') 8, KSEC2(8)
      write(*,'(10x,"KSEC2date",I7,"-",I2,"-",I2,"_",I2)') KSEC2(10),KSEC2(11),KSEC2(12),KSEC2(13)
      endif
      else
        write(*,'(/,8x,"Error MM5toGRIB: Not pressure levels",/)')
        stop
      endif
   end subroutine header_grib


subroutine read_namelist_output(ltogrib,nametogrib,lmore,sksec)

 implicit none

 character(len=20) :: listfile='namelist.output'
 character(len=200) :: listline
 character(len=20) :: parametre, parametreb
 Character(len=1) :: parametrec

 integer :: iunit=15
 integer :: ierr

 character (len= 9), dimension(200) :: nametogrib
 logical, dimension(200) :: ltogrib
 logical :: lmore

 integer :: in
! integer :: inn_end
! parameter (inn_end = 200)

 type sksecone
         integer :: six
         integer :: seven
         integer :: eight 
         integer :: nine 
         real    :: factor
 end type sksecone

 type(sksecone),dimension(200) :: sksec
 integer :: ic, ici, icf

 nametogrib='         '

 open(iunit, file=listfile, err=599, form='formatted', status='old',&
      action='read')

      if (lmore) write(*,'(/,8x,"Subroutine read_namelist_output: Reading namelist.output file ...",/)')
      read(iunit, '(A100)', iostat=ierr, end=500) listline
      if (lmore) write(*,'(A100)') listline
      read(iunit, '(A100)', iostat=ierr, end=500) listline
      if (lmore) write(*,'(A100)') listline
      read(iunit, '(A100)', iostat=ierr, end=500) listline
      if (listline(1:8).ne.'&record1') then
         parametre='&record1'
         goto 596
      endif
      if (lmore) write(*,'(A100)') listline

      in=1
      do while (ierr == 0)
           read(iunit, 510, iostat=ierr, end=500,err=594) parametre, nametogrib(in), parametreb, ltogrib(in), listline
           if ((parametre(1:8).ne.' name: #').and.(parametreb(1:12).ne.'#  togrib: .')) goto 597
           if ((ltogrib(in)).and.(lmore)) print*, parametre(1:8),nametogrib(in),parametreb(1:12),ltogrib(in),listline
           ic=0
           ici=0
           icf=0
!           print*, listline

           read(listline(ic:ic),'(A1)',err=594) parametrec 
!           print*, parametrec
           do while (parametrec.ne.'&')
             ic=ic+1
             if (ic.eq.200) goto 595
             read(listline(ic:ic),'(A1)',err=594) parametrec
!             print*, parametrec
           enddo
           ic=ic+1
           ici=ic
           read(listline(ic:ic),'(A1)',err=594) parametrec
           do while (parametrec.ne.'&')
             ic=ic+1
             if (ic.eq.200) goto 595
             read(listline(ic:ic),'(A1)',err=594) parametrec
!             print*, parametrec
           enddo
           ic=ic+1
           icf=ic-2
!           print*, listline(ici:icf)
           read(listline(ici:icf),'(I100)',err=594) sksec(in)%six
           ici=ic
           read(listline(ic:ic),'(A1)',err=594) parametrec
           do while (parametrec.ne.'&')
             ic=ic+1
             if (ic.eq.200) goto 595
             read(listline(ic:ic),'(A1)',err=594) parametrec
!             print*, parametrec
           enddo
           ic=ic+1
           icf=ic-2
!           print*, listline(ici:icf)
           read(listline(ici:icf),'(I100)',err=594) sksec(in)%seven
           ici=ic
           read(listline(ic:ic),'(A1)',err=594) parametrec
           do while (parametrec.ne.'&')
             ic=ic+1
             if (ic.eq.200) goto 595
             read(listline(ic:ic),'(A1)',err=594) parametrec
!             print*, parametrec
           enddo
           ic=ic+1
           icf=ic-2
!           print*, listline(ici:icf)
           read(listline(ici:icf),'(I100)',err=594) sksec(in)%eight
           ici=ic
           read(listline(ic:ic),'(A1)',err=594) parametrec
           do while (parametrec.ne.'&')
             ic=ic+1
             if (ic.eq.200) goto 595
             read(listline(ic:ic),'(A1)',err=594) parametrec
!             print*, parametrec
           enddo
           ic=ic+1
           icf=ic-2
!           print*, listline(ici:icf)
           read(listline(ici:icf),'(I100)',err=594) sksec(in)%nine
           ici=ic
           read(listline(ici:200),'(F50.40)',err=594) sksec(in)%factor
!           print*, sksec(in)

!           print*, in
!           if (in.eq.4) stop 
           in=in+1
!           stop
      enddo
      stop
      if (lmore) print*, ' '

510 format(A8,A9,A12,L1,A100)
      goto 500

594   write(*,'(/,8x, "Subroutine read_namelist_output: Error reading in ", A,": ",/,10x,A8,A9,A12,L1,A200,/)')&
               trim(listfile), trim(parametre), trim(nametogrib(in)), trim(parametreb), ltogrib(in), listline
      stop 594

595   write(*,'(/,8x, "Subroutine read_namelist_output: Error reading in ", A,": ",/,10x,A8,A9,A12,L1,A200,/)')&
               trim(listfile), trim(parametre), trim(nametogrib(in)), trim(parametreb), ltogrib(in), listline
      stop 595

596   write(*,'(/,8x, "Subroutine read_namelist_output: Error reading in", A, 1x, "the parameter ", A)')&
               trim(listfile), trim(parametre)
      stop 596
597   write(*,'(/,8x, "Subroutine read_namelist_output: Error reading in", A, 1x, ": ",A8,A9,A12,L1,".")')&
               trim(listfile), trim(parametre), trim(nametogrib(in)), trim(parametreb), ltogrib(in)
      stop 597
500   close(iunit, err=598)


      goto 600
598   write(*,'(/,8x, "Subroutine read_namelist_output: Error closing the ", A, 1x, "file",/)')&
               trim(listfile)
      stop 598
599   write(*,'(/,8x, "Subroutine read_namelist_output: Cannot open the ", A, 1x, "file",/)')&
               trim(listfile)
      stop 599
600   continue

!      return

end subroutine read_namelist_output


