subroutine output(hdate, nlvl, maxlvl, plvl, interval, iflag, out_format, prefix, debug_level)
!                                                                             !
!*****************************************************************************!
!  Write output to a file.
!                                                                             !
!    hdate:
!    nlvl:  
!    maxlvl:
!    plvl:
!    interval:
!    iflag:  1 = output for ingest into rrpr ; 2 = output for hinterp
!    out_format: requested output format (WPS, SI, or MM5)
!                                                                             !
!*****************************************************************************!

  use table
  use gridinfo
  use storage_module
  use filelist
  use module_debug
  use stringutil

  implicit none

  character(LEN=19) :: hdate
  character(LEN=24) :: hdate_output
  character(LEN=3)  :: out_format
  character(LEN=256)  :: prefix
  integer :: iunit = 13

  real, pointer, dimension(:,:) :: scr2d

  integer :: maxlvl
  integer nlvl, debug_level
  real , dimension(maxlvl) :: plvl
  character (LEN=9) :: field
  real :: level
  integer :: sunit = 14
  integer :: interval
  integer :: iflag
! Local Miscellaneous
  integer :: k, n, m, ilev
  integer :: ii, jj
  real :: maxv, minv
  real :: xplv
  real :: xfcst = 0.
  character (LEN=25) :: units
  character (LEN=46) :: Desc
  logical lopen

! DATELEN:  length of date strings to use for our output file names.
  integer :: datelen

! Decide the length of date strings to use for output file names.  
! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
  if (mod(interval,3600) == 0) then
     datelen = 13
  elseif (mod(interval,60) == 0) then
     datelen = 16
  else
     datelen = 19
  endif
 
  call get_plvls(plvl, maxlvl, nlvl)

  if ( debug_level .ge. 0 ) then
  write(*,119) hdate(1:10), hdate(12:19)
119 format(/,79('#'),//,'Inventory for date = ', A10,1x,A8,/)

  write(*,advance='NO', fmt='("PRES", 2x)')
  WRTLOOP : do n = 1, maxvar
     do k = 1, n-1
        if (namvar(k).eq.namvar(n)) cycle WRTLOOP
     enddo
     write(*,advance='NO', fmt='(1x,A8)') namvar(n)
  enddo WRTLOOP
  write(*,advance='YES', fmt='(1x)')

  write(*,FMT='(79("-"))')
  end if
  KLOOP : do k = 1, nlvl
     if ((iflag.eq.2).and.(plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
        cycle KLOOP
     endif
     ilev = nint(plvl(k))
     if ( debug_level .ge. 0 ) then
     write(*, advance='NO', FMT='(F6.1)') plvl(k)/100.
     end if
     MLOOP : do m = 1, maxvar
        do n = 1, m-1
           if (namvar(m).eq.namvar(n)) cycle MLOOP
        enddo
        if ( debug_level .ge. 0 ) then
        if (is_there(ilev,namvar(m))) then
           write(*, advance='NO', FMT='("  X      ")')
        else
	   if ( plvl(k).gt.200000 ) then
             write(*, advance='NO', FMT='("  O      ")')
	   else
             write(*, advance='NO', FMT='("         ")')
	   endif
        endif
        endif
     enddo MLOOP
     if ( debug_level .ge. 0 ) then
     write(*,advance='YES', fmt='(1x)')
     endif
  enddo KLOOP
  if ( debug_level .ge. 0 ) then
  write(*,FMT='(79("-"))')
  endif
  
  if (iflag.eq.1) then
     if (nfiles.eq.0) then
        open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
             position='REWIND')
        nfiles = nfiles + 1
        filedates(nfiles)(1:datelen) = hdate(1:datelen)
     else
        DOFILES : do k = 1, nfiles
           if (hdate(1:datelen).eq.filedates(k)(1:datelen)) then
              open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted',&
                   position='APPEND')
           endif
        enddo DOFILES
        inquire (iunit, OPENED=LOPEN)
        if (.not. LOPEN) then
           open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
                position='REWIND')
           nfiles = nfiles + 1
           filedates(nfiles)(1:datelen) = hdate(1:datelen)
        endif
     endif
  else if (iflag.eq.2) then
     open(iunit, file=trim(prefix)//':'//HDATE(1:datelen), form='unformatted', &
          position='REWIND')
  endif

!MGD  if ( debug_level .gt. 100 ) then
!MGD     write(6,*) 'begin nloop'
!MGD  end if
  NLOOP : do n = 1, nlvl

!MGD  if ( debug_level .gt. 100 ) then
!MGD     write(6,*) 'begin outloop'
!MGD  end if
     OUTLOOP : do m = 1, maxvar
        field = namvar(m)
        do k = 1, m-1
           if (field.eq.namvar(k)) cycle OUTLOOP
        enddo
        level = plvl(n)
        if ((iflag.eq.2).and.(level.gt.200100) .and. (level.lt.200200)) then
           cycle NLOOP
        endif
        ilev = nint(level)
        desc = ddesc(m)
        if (iflag.eq.2) then
           if (desc.eq.' ') cycle OUTLOOP
        endif
        units = dunits(m)
        if ((iflag.eq.1).or.(iflag.eq.2.and.desc(1:1).ne.' ')) then
          if (is_there(ilev,field)) then 
            call get_dims(ilev, field)

!MGD            if ( debug_level .gt. 100 ) then
!MGD               write(6,*) 'call refr_storage'
!MGD            end if
            call refr_storage(ilev, field, scr2d, map%nx, map%ny)

!MGD            if ( debug_level .gt. 100 ) then
!MGD               write(6,*) 'back from refr'
!MGD               write(6,*) 'out_format = ',out_format
!MGD            end if

	    if (out_format(1:2) .eq. 'SI') then
!MGD              if ( debug_level .gt. 100 ) then
!MGD                 write(6,*) 'writing in SI format'
!MGD              end if
              write(iunit) 4
              hdate_output = hdate
              write (iunit) hdate_output, xfcst, map%source, field, units, &
                   Desc, level, map%nx, map%ny, map%igrid
              if (map%igrid.eq.3) then ! lamcon
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
                      map%lov, map%truelat1, map%truelat2
              elseif (map%igrid.eq.5) then ! Polar Stereographic
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
                      map%lov, map%truelat1
              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
              elseif (map%igrid.eq.1)then ! Mercator
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
                      map%truelat1
              else
                 call mprintf(.true.,ERROR, &
                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
              endif
              write (iunit) scr2d
	    else if (out_format(1:2) .eq. 'WP') then   
                call mprintf(.true.,DEBUG, &
         "writing in WPS format  iunit = %i, map%%igrid = %i",i1=iunit,i2=map%igrid)
              write(iunit) 5
              hdate_output = hdate
              write (iunit) hdate_output, xfcst, map%source, field, units, &
                   Desc, level, map%nx, map%ny, map%igrid
              if (map%igrid.eq.3) then ! lamcon
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
                      map%lov, map%truelat1, map%truelat2, map%r_earth
              elseif (map%igrid.eq.5) then ! Polar Stereographic
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
                      map%lov, map%truelat1, map%r_earth
              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
		      map%r_earth
              elseif (map%igrid.eq.1)then ! Mercator
                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
                      map%truelat1, map%r_earth
              else
                 call mprintf(.true.,ERROR, &
                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
              endif
	      write (iunit) map%grid_wind
              write (iunit) scr2d
	    else if (out_format(1:2) .eq. 'MM') then
!MGD              if ( debug_level .gt. 100 ) then
!MGD	         write(6,*) 'writing in MM5 format'
!MGD              end if
              write(iunit) 3
              hdate_output = hdate
              write (iunit) hdate_output, xfcst, field, units, Desc, level,&
                   map%nx, map%ny, map%igrid
              if (map%igrid.eq.3) then ! lamcon
                 write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
                      map%truelat1, map%truelat2
              elseif (map%igrid.eq.5) then ! Polar Stereographic
                 write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
                      map%truelat1
              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
                 write (iunit) map%lat1, map%lon1, map%dy, map%dx
              elseif (map%igrid.eq.1)then ! Mercator
                 write (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
              else
                 call mprintf(.true.,ERROR, &
                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
              endif
              write (iunit) scr2d
	    endif
              if ( debug_level .gt. 100 ) then
	        call mprintf(.true.,DEBUG, &
	        "hdate = %s,  xfcst = %f ",s1=hdate_output,f1=xfcst)
	        call mprintf(.true.,DEBUG, &
           "map%%source = %s, field = %s, units = %s",s1=map%source,s2=field,s3=units)
	        call mprintf(.true.,DEBUG, &
	            "Desc = %s, level = %f",s1=Desc,f1=level)
	        call mprintf(.true.,DEBUG, &
	            "map%%nx = %i, map%%ny = %i",i1=map%nx,i2=map%ny)
              else if ( debug_level .gt. 0 ) then
	        call mprintf(.true.,STDOUT, &
	      " field = %s, level = %f",s1=field,f1=level)
	        call mprintf(.true.,LOGFILE, &
	      " field = %s, level = %f",s1=field,f1=level)
              end if
              if ( debug_level .gt. 100 ) then
	      maxv = -99999.
	      minv = 999999.
	      do jj = 1, map%ny
	      do ii = 1, map%nx
	        if (scr2d(ii,jj) .gt. maxv) maxv = scr2d(ii,jj)
	        if (scr2d(ii,jj) .lt. minv) minv = scr2d(ii,jj)
	      enddo
	      enddo
	      call mprintf(.true.,DEBUG, &
	         "max value = %f , min value = %f",f1=maxv,f2=minv)
              end if

              nullify(scr2d)

           endif
        endif
     enddo OUTLOOP
  enddo NLOOP

  close(iunit)

end subroutine output

