program merge 
c
c**   This program combines the poe probabilities for multiple storms
c**     into a single probability array to create GrADS and GRIB files
c
      character*1  poe_num
      character*10 cur_ymdh
c
c**   Read in the current DTG and the total number of poe processes 
c**    from the command line
c
      call getarg ( 1, cur_ymdh )
      call getarg ( 2, poe_num )
c
c**   Generate the final probabilities from the poe output files
c
      call final_prblty ( poe_num, cur_ymdh )
c
c      if ( poe_num .eq. 'X' ) stop
c     & ' ***** There are no probabilities to process! *****'
c
c**   Prepare the GrADS plot file
c
      call plotter ( cur_ymdh )
c
c**   Prepare the GRIB1 plot file
c
      call gribout ( cur_ymdh )
c
c**   Prepare the GRIB2 plot file
c
      call grib2out ( cur_ymdh )
c
      stop ' *** The poe final probability program has finished! ***'
c
      end
c
c*************************************************************************
      subroutine final_prblty ( poe_num, cur_ymdh )
c
c**     This program computes the Demaria probabilities
c
c**   Parameter statement for 1.0 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=15,imax=110,ioffset=-14,iarray=96,deltax=1.0)
cc      parameter (jmin= 1,jmax= 60,joffset=  0,jarray=60,deltay=1.0)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=30,imax=220,ioffset=-29,iarray=191,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
cc      parameter (ix=22,jx=12)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 100.0E
c
cc      parameter (imin=30,imax=520,ioffset=-29,iarray=491,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
     
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 1.0W to 100.0E
c
      parameter (imin=2,imax=520,ioffset=-1,iarray=519,deltax=0.5)
      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.1 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=150,imax=1100,ioffset=-149,iarray=951,deltax=0.1)
cc      parameter (jmin= 10,jmax= 600,joffset=  -9,jarray=591,deltay=0.1)
c
cc      parameter ( numits=10 )
      parameter ( numits=20 )
      parameter ( numitsh=10 )
c
      common /final/ parray( iarray, jarray, 0:numits, 3 ),
     &               parrayi( iarray, jarray, 0:numits, 3 ),
     &               parrayh( iarray, jarray, 0:numitsh, 3 )
c
      logical file_status
c
      character*1  poe_num, poe_char
      character*10 cur_ymdh
      character*20 filename
      real left,right,top,bottom
c
      data file_status / .false. /
c
c**   set the final probability array to the default value
c
      do n = 1, 3
         do k = 0, numits
            do j = 1, jarray
               do i = 1, iarray 
                  parray( i, j, k, n ) = 0.0
                  parrayi( i, j, k, n ) = 0.0
               enddo
            enddo
         enddo
      enddo
c
c**   Open the process probability files, read in the values and
c**     accumulate them in the final probability array
c
      read ( poe_num, '(i1)' ) num_poe
c
      do num = 1, num_poe
c
         bottom = 9999 
         top = -9999
         left = -9999
         right = 9999
	 
         write ( poe_char, '(i1)' ) num
         filename = 'prblty_'//poe_char//'.dat'
c
         open ( 21, file=filename, status='old', form='unformatted',
     &        err=10 )
c
         do n = 1, 3
            do k = 0, numits
               do j = jmax, jmin, -1
                  do i = imax, imin, -1
c 
                     read (21,iostat=ios,err=1010) temp_prob, temp_probi
c
c**   Probabilities must be added in a special way:
c**       P(ab) = P(a) + P(b) - P(a)*P(b)/100
c
          prblty = parray( i + ioffset, j + joffset, k, n )
          prblty = prblty + temp_prob - prblty*temp_prob/100.0
          parray( i + ioffset, j + joffset, k, n ) = prblty 
c
          prblty = parrayi( i + ioffset, j + joffset, k, n )
          prblty = prblty + temp_probi - prblty*temp_probi/100.0
          parrayi( i + ioffset, j + joffset, k, n ) = prblty 

c**   Find corners of the grid, for each storm	  
	  if(temp_prob .ge. 1 .and. k .eq. numits) then	  	
		if(i .lt. right) right = i
		if(i .gt. left) left = i
		
		if(j .lt. bottom) bottom = j
		if(j .gt. top) top = j
	  endif
c
                  enddo
               enddo
            enddo
         enddo
c
         file_status = .true.

c** calculate corners	 
	 right = 0 - (right * .5)
	 left = 0 - (left * .5)
	 top = top * .5
	 bottom = bottom * .5

c**  write out corners
         filename = cur_ymdh//'.corners'
	 open ( 22, file=filename, status='old', form='formatted',
     &        position = 'APPEND', err=10 )
     
         write(22,*) bottom,';',left,';',top,';',right
c
	 close(22)
c
c
 10      continue
c
         close ( 21 )
c
      enddo
c


c     Compute 12 hourly incremental from 6 hourly cumulative and 6 hourly incremental
c     using the formula
c         i12[t] = i6[t-6] + (c6[t] - c6[t-6])
c     where t is the time period (12-120 by 12), i12 is the twelve hourly incremental 
c     probability, i6 is the six hourly incremental probability, and c6 is the 6 hourly
c     cumulative probability
      do n = 1,3
        do k = 0, numitsh
	  
	  if (k==0) then
	    ki = 0
            kj = 0
	  else
	      ki = k * 2
              kj = ki - 1 
	  endif
	      
          do j = jmax, jmin, -1
            do i = imax, imin, -1      	
	      
	      parrayh( i + ioffset, j + joffset, k, n) =
     &          parrayi( i + ioffset, j + joffset, kj, n) + (
     &          parray( i + ioffset, j + joffset, ki, n) - 
     &          parray( i + ioffset, j + joffset, kj, n)     
     &          )
     
            enddo
	  enddo
	enddo
      enddo



      if ( .not. file_status ) poe_num = 'X'
c
      return
c
 1010 print *, '*** ERROR - reading file = ', filename, ' ios = ', ios
      stop
c
      end
c
c********1*********2*********3*********4*********5*********6*********7**
      subroutine plotter ( ymdh )
C
C**   This program prepares the grads plot file for the radii array.
c
c**   Parameter statement for 1.0 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=15,imax=110,ioffset=-14,iarray=96,deltax=1.0)
cc      parameter (jmin= 1,jmax= 60,joffset=  0,jarray=60,deltay=1.0)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=30,imax=220,ioffset=-29,iarray=191,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 100.0E
c
cc      parameter (imin=30,imax=520,ioffset=-29,iarray=491,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
      
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 1.0W to 100.0E
c
      parameter (imin=2,imax=520,ioffset=-1,iarray=519,deltax=0.5)
      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.1 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=150,imax=1100,ioffset=-149,iarray=951,deltax=0.1)
cc      parameter (jmin= 10,jmax= 600,joffset=  -9,jarray=591,deltay=0.1)
c
cc      parameter ( numits=10 )
      parameter ( numits=20 )
      parameter ( numitsh=10 )
c
c      common /final/ parray( iarray, jarray, 0:numits, 3 ),
c     &     parrayi( iarray, jarray, 0:numits, 3 )
      common /final/ parray( iarray, jarray, 0:numits, 3 ),
     &               parrayi( iarray, jarray, 0:numits, 3 ),
     &               parrayh( iarray, jarray, 0:numitsh, 3 )
c
      integer end_plot
c
      character*3  month(12)
      character*10 ymdh
      character*12 time_line
c
      data month / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN',
     &             'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' /
c
c**   Create the time line
c
      read ( ymdh(5:6), '( i2)' ) nmonth
      time_line = ymdh(9:10)//'Z'//ymdh(7:8)//month(nmonth)//ymdh(1:4)
c
      xstart = - imax*deltax
      ystart =   jmin*deltay
c
c**   Open, write and close the wind grads control file
c
      open ( 31, file='prob.ctl', status='unknown')
c
      write ( 31, '( ''dset   prob.gr'' )' )
      write ( 31, '( ''title  Tropical Cyclone Probability Field'' )' )
      write ( 31, '( ''undef  -999.0'' )' )
      write ( 31, '( ''xdef   '', i4, '' linear'', 2f7.1 ) ') 
     &               iarray, xstart, deltax
      write ( 31, '( ''ydef   '', i4, '' linear'', 2f7.1 ) ') 
     &               jarray, ystart, deltay
      write ( 31, '( ''zdef      1 levels   1   1'' )' )
cc      write ( 31, '( ''tdef     '', i2,'' linear '', a12,'' 12hr'' )') 
cc     &               numits + 1, time_line
      write ( 31, '( ''tdef     '', i2,'' linear '', a12,'' 6hr'' )') 
     &               numits + 1, time_line
      write ( 31, '( ''*'' ) ')
      write ( 31, '( ''vars      6'' ) ')
      write ( 31, '( ''p_34      1 99  34-knot tprobability field'' ) ')
      write ( 31, '( ''p_50      1 99  50-knot tprobability field'' ) ')
      write ( 31, '( ''p_64      1 99  64-knot tprobability field'' ) ')
      write ( 31, '( ''pi_34     1 99  34-knot iprobability field'' ) ')
      write ( 31, '( ''pi_50     1 99  50-knot iprobability field'' ) ')
      write ( 31, '( ''pi_64     1 99  64-knot iprobability field'' ) ')
      write ( 31, '( ''endvars'' ) ')
c
      close ( 31 )
c
c**   Open, write and close the wind grads file
c
      open ( 32, file='prob.gr', status='unknown', 
     &     access='direct', form='unformatted', recl=iarray*jarray*4 )
c
      ict = 0
c
      do n = 0, numits
c
c**   Write out the 34-knot probability field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &       ((parray( i, j, n, 1 ), i = iarray, 1, -1 ), j = 1, jarray)
c
c**   Write out the 50-knot probability field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &       ((parray( i, j, n, 2 ), i = iarray, 1, -1 ), j = 1, jarray)
c
c**   Write out the 64-knot probability field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &       ((parray( i, j, n, 3 ), i = iarray, 1, -1 ), j = 1, jarray)
c
c**   Write out the 34-knot probability incremental field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &      ((parrayi( i, j, n, 1 ), i = iarray, 1, -1 ), j = 1, jarray)
c
c**   Write out the 50-knot probability incremental field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &      ((parrayi( i, j, n, 2 ), i = iarray, 1, -1 ), j = 1, jarray)
c
c**   Write out the 64-knot probability incremental field
c
          ict = ict + 1
c
          write ( 32, rec=ict ) 
     &      ((parrayi( i, j, n, 3 ), i = iarray, 1, -1 ), j = 1, jarray)
c
      enddo
c
      close ( 32 )
c
      return
      end
c********1*********2*********3*********4*********5*********6*********7**
      subroutine gribout ( ymdh )
c
c**   The purpose of this subroutine is to create an output grib message
c**    for each forecast hour.  The subroutine creates the PDS and GDS 
c**    sections of the GRIB header, then calls the NCEP putgb routine to
c**    pack the GRIB message.
c
c**   Parameter statement for 1.0 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=15,imax=110,ioffset=-14,iarray=96,deltax=1.0)
cc      parameter (jmin= 1,jmax= 60,joffset=  0,jarray=60,deltay=1.0)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=30,imax=220,ioffset=-29,iarray=191,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 100.0E
c
cc      parameter (imin=30,imax=520,ioffset=-29,iarray=491,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)

c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 1.0W to 100.0E
c
      parameter (imin=2,imax=520,ioffset=-1,iarray=519,deltax=0.5)
      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)

c
c**   Parameter statement for 0.1 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=150,imax=1100,ioffset=-149,iarray=951,deltax=0.1)
cc      parameter (jmin= 10,jmax= 600,joffset=  -9,jarray=591,deltay=0.1)
c
      integer numits, deltat
c
cc      parameter ( numits = 10, deltat = 12 )
      parameter ( numits = 20, deltat = 6 )
      parameter ( numitsh = 10  )
c
C      common /final/ parray( iarray, jarray, 0:numits, 3 ),
C     &               parrayi( iarray, jarray, 0:numits, 3 )
	common /final/ parray( iarray, jarray, 0:numits, 3 ),
     &               parrayi( iarray, jarray, 0:numits, 3 ),
     &               parrayh( iarray, jarray, 0:numitsh, 3 )
c
      real       temparr(iarray,jarray)
      integer    itype1( 3 ), itype2 ( 3 ), kpds( 25 ), kgds( 22 )
      integer    end_plot, ngribunit
      logical(1) lb(iarray,jarray)
c
      character*10 ymdh
c
      data itype1 / 193, 194, 195 /
      data itype2 / 198, 199, 200 /
      data ngribunit / 61 /

cc      lb = .TRUE.

      read (ymdh(1:4), '(i4)') iyyyy
      read (ymdh(1:2), '(i2)') icc
      read (ymdh(3:4), '(i2)') iyy
      read (ymdh(5:6), '(i2)') imm
      read (ymdh(7:8), '(i2)') idd
      read (ymdh(9:10),'(i2)') ihh

      icent = icc + 1

      print *,'iyyyy= ',iyyyy,'  imm= ',imm,' idd= ',idd,' ihh= ',ihh
      print *,'ymdh= ',ymdh

c     A GRIB file is made up of 1 or more individual GRIB records.
c     Each individual GRIB record is made up of 4 sections that 
c     you need to worry about: 
c
c     (1) Product Definition Section (PDS): This contains 
c         integer values that define values to describe the type
c         of data that will be stored in this record, where it 
c         came from (TPC), what the starting date is, what the 
c         forecast hour is, what vertical level the data are on,
c         what precision you would like preserved for the data
c         values, etc....
c
c     (2) Grid Description Section (GDS): This contains integer
c         values that describe the map projection of the grid, 
c         the grid spacing, the beginning and end points of the 
c         grid, and flags that indicate if your data values go
c         from north to south, east to west, or vice versa.
c         Technically, if you are using a standard grid that is
c         known and identified in the NCEP GRIB guide (such as 
c         one of the Eta or AVN grids), you could just enter the
c         ID number for that grid in kpds(3) in the PDS and not 
c         even have to create a GDS.  However, it is NCEP's 
c         practice to include a GDS for all GRIB records that 
c         they produce; it makes it easier for the end user to
c         know what you're doing.
c
c     (3) Bit Map Section (BMS): This is a small section made 
c         up of one-bit logical values.  There is one logical 
c         value for each point on your grid, and it indicates 
c         whether or not valid data exists at each gridpoint.  
c         It is most often used for applications such as when 
c         packing sea surface temperature data, and you want 
c         to mask out land data points.  You would NOT want 
c         to use this, for example, on a rainfall map to 
c         indicate areas that don't have any accumulated 
c         rainfall, or on your wind maps to indicate areas 
c         not included within the 34-knot wind areas.  It is
c         optional to put this section in.  I usually always 
c         include it and turn all values to .TRUE. by default.
c
c     (4) Binary Data Section (BDS): The data values are 
c         packed as a series of binary digits.  One critical 
c         step in this process that the GRIB software takes 
c         care of for you is that it scans the data, finds 
c         the minimum value, and subtracts that value from 
c         every value in your array.  This has 2 benefits:
c         it removes negative numbers and any machine-specific
c         problems that may imply, and for data values that 
c         are significantly far from 0 (such as pressure in 
c         Pascals), it can significantly compress the space 
c         needed to store the data.  The reference value is 
c         stored and automatically added to all the values 
c         by the software when being read by a user.


c     Set the PDS section of the GRIB header.  Note that kpds(14)
c     is the forecast hour and is set in the loop below.  Also, 
c     kpds(22) is the "units decimal scale factor".  Essentially,
c     it tells the GRIB software how many digits of precision you
c     want to keep after the decimal point (positive integer, for
c     something like vorticity values) or before the decimal 
c     point (negative integer, for something like storing pressure
c     in Pascals (eg, for 102042 Pa, entering a -1 for kpds(22) 
c     would mean that when you unpacked, or read, this record, 
c     you would read a value of 10204 Pa (1020.4 mb) for this 
c     record)).  For your wind values, I set it to 0, since you
c     only care about whole number wind values.

      kpds(1)  =      52  ! = forecast center is TPC/NHC
      kpds(2)  =      30  ! = "Forecaster generated field"
      kpds(3)  =     255  ! = grid is user-defined in GDS
      kpds(4)  =     192  ! = GDS & BMS both included in header
cc      kpds(5)  =     194  ! = parameter is probability of snow
      kpds(6)  =       1  ! = "type of level" code for surface
      kpds(7)  =       0  ! = actual height level is surface level
      kpds(8)  =     iyy  ! = initial yyyy of forecast
      kpds(9)  =     imm  ! = initial   mm of forecast
      kpds(10) =     idd  ! = initial   dd of forecast
      kpds(11) =     ihh  ! = initial   hh of forecast
      kpds(12) =       0  ! = initial minutes value of forecast
      kpds(13) =       1  ! = value of kpds(14) is in hours
C      kpds(14) =       0  ! Fcst hour defined in loop below
      kpds(15) =       0  ! = 2nd fcst hour, N/A for you....
      kpds(16) =      10  ! = fcst valid at ref time + kpds(14)
      kpds(17) =       0  ! = # included in average, N/A for you
      kpds(18) =       1  ! = version # of GRIB
      kpds(19) =       2  ! = version # of GRIB parameter table
      kpds(20) =       0  ! = # missing from avg, N/A for you
      kpds(21) =   icent  ! = century of reference time of data
      kpds(22) =       1  ! = units decimal scale factor
cc      kpds(22) =       0  ! = units decimal scale factor
      kpds(23) =       0  ! = forecast sub-center, N/A
      kpds(24) =       0  ! = for ensemble products use only
      kpds(25) =       0  ! = not used

c     Set the GDS section of the GRIB header.  kpds(11) is the 
c     scanning mode flag, which is an integer number that indicates
c     how the data scan, i.e., in the +i,-j direction, or in the
c     -i,+j direction, etc....  A value of 0 indicates that the
c     scanning is in the +i,-j direction, implying that the data
c     values start in the upper left corner of the grid and end 
c     in the lower right corner. 

      nxstart = (nint (360.0 - (float(iarray) * deltax) + 
     &                        (float(ioffset) * deltax))) * 1000
      nyend   = (nint (float(jmax) * deltay)) * 1000

      nxend   = nxstart + (nint (float(iarray-1) * deltax)) * 1000
      nystart = nyend   - (nint (float(jarray-1) * deltay)) * 1000

      ndeltax = nint (deltax * 1000.0)
      ndeltay = nint (deltay * 1000.0)

      kgds(1)  =       0  ! = Data representation type (0 = lat/lon)
      kgds(2)  =  iarray  ! = # pts in x-direction
      kgds(3)  =  jarray  ! = # pts in y-direction
      kgds(4)  = nystart  ! = Latitude of first point (*1000)
      kgds(5)  = nxstart  ! = Longitude of first point (*1000)
      kgds(6)  =     128  ! = resolution & components flag
      kgds(7)  =   nyend  ! = Latitude of last point (*1000)
      kgds(8)  =   nxend  ! = Longitude of last point (*1000)
      kgds(9)  = ndeltax  ! = x-increment in degrees * 1000
      kgds(10) = ndeltay  ! = y-increment in degrees * 1000
      kgds(11) =      64  ! = scanning mode flag
      kgds(12) =       0  ! 12-19 N/A for you.....
      kgds(13) =       0      
      kgds(14) =       0
      kgds(15) =       0      
      kgds(16) =       0
      kgds(17) =       0      
      kgds(18) =       0
      kgds(19) =       0      
      kgds(20) =     255  ! = N/A for you, leave as 255....
c
      kf = iarray * jarray
c
c**   Now run through for each forecast hour and call the NCEP
c**   w3lib routine "putgb", which packs the data into 1 output
c**   GRIB record and writes it to unit "ngribunit".
c
      call open_grib_files ( ngribunit, iret )
c
      if (iret /= 0) then
         print '(/,a50,2i4,/)', '!!! ERROR: opening unit, rc= ', 
     &        ngribunit, iret
         stop 99
      endif
c
cc      open  ( 31, file='grib1_input_array', status='unknown' )
c
      do n = 0, numits
c
         kpds(14) = n*deltat

c
c**   Write out the 34-knot, 50-knot and 64-knot cummulative probabilities
c
         do m = 1, 3
c
            do j = 1, jarray
               do i = 1, iarray
c
c**   Flip the probability array in the x-direction
c
                  iflip = iarray - i + 1
                  temparr( iflip, j ) = parray( i, j, n, m )
c
c**   Set the bit map for valid data
c
                  lb( iflip, j ) = .true.
c  		
               enddo
            enddo
c
cc            write ( 31, '( 2i6 )' ) kpds(14), m
cc            write ( 31, '(119f5.1)' ) 
cc     &           ((temparr( i, j ), j = 1, jarray), i = 1, iarray)
c
            kpds( 5 ) = itype1( m ) 
c
            call putgb ( ngribunit, kf, kpds, kgds, lb, temparr, iret )
c
            if ( iret == 0 ) then
               print *, ' '
               print *, '+++ IRET = 0 after call to putgb'
               print *, ' '
            else
               print *, ' '
               print *, '!!!IRET NE 0 AFTER CALL TO PUTGB!!!'
               print *, '!!!IRET = ', iret
               print *, ' '
            endif
c
         enddo
c
c**   Write out the 34-knot, 50-knot and 64-knot incremental probabilities
c
         do m = 1, 3
c
            do j = 1, jarray
               do i = 1, iarray
c
c**   Flip the probability array in the x-direction
c
                  iflip = iarray - i + 1
                  temparr( iflip, j ) = parrayi( i, j, n, m )
c
c**   Set the bit map for valid data
c
                  lb( iflip, j ) = .true.
c  		
               enddo
            enddo
c
            kpds( 5 ) = itype2( m ) 
c
            call putgb ( ngribunit, kf, kpds, kgds, lb, temparr, iret )
c
            if ( iret == 0 ) then
               print *, ' '
               print *, '+++ IRET = 0 after call to putgb'
               print *, ' '
            else
               print *, ' '
               print *, '!!!IRET NE 0 AFTER CALL TO PUTGB!!!'
               print *, '!!!IRET = ', iret
               print *, ' '
            endif
c
         enddo
c
      enddo

cc      close ( 31 )
c
      call baclose ( ngribunit, iret1 )
c
      print *, 'baclose return code for grib1out: iret1= ', iret1
c
      return
      end
c
c********1*********2*********3*********4*********5*********6*********7**
      subroutine open_grib_files ( unit1, iret )
C
C     ABSTRACT: This subroutine must be called before any attempt is
C     made to read from the input GRIB files.  The GRIB and/or index 
C     files are opened with a call to baopenr/baopenw.  These calls to 
C     the baopen routines were not needed in the cray version of this 
C     program (the files could be opened with a simple Cray assign 
C     statement), but the GRIB-reading utilities on the SP do require 
C     calls to these routines (it has something to do with the GRIB 
C     I/O being done in C on the SP, and the C I/O package needs an 
C     explicit open statement).  NOTE: Input GRIB files are opened 
C     with a call to baopenr, while output GRIB files are opened with
C     a call to baopenw.
C
C     INPUT:
C     unit1   The Fortran unit number for the GRIB data
c             and/or GRIB index files
C
C     OUTPUT:
C     iret    The return code from this subroutine

      character fname*7
      integer   unit1
c
c**   Open the output grib file....
c
      fname(1:5) = "fort."
      write ( fname(6:7), '(I2)' ) unit1
c
      call baopenw ( unit1, fname, iret1 )
c      
      print *, 'baopenw (unit 51): iret1= ', iret1
c      
      if ( iret1 /= 0 ) then
        print *, ' '
        print *, '!!! ERROR in sub open_grib_files opening grib'
        print *, '!!! or grib index file.  baopen return codes:'
        print *, '!!! unit 51 baopenw return code = iret1 = ', iret1
        iret = 93
        return
      endif
        iret = 0
c        
      return
      end
c********1*********2*********3*********4*********5*********6*********7**
      subroutine grib2out ( ymdh )
c
c**   The purpose of this subroutine is to create an output grib message
c**    for each forecast hour.  The subroutine creates the PDS and GDS 
c**    sections of the GRIB header, then calls the NCEP putgb routine to
c**    pack the GRIB message.
c
c**   Parameter statement for 1.0 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=15,imax=110,ioffset=-14,iarray=96,deltax=1.0)
cc      parameter (jmin= 1,jmax= 60,joffset=  0,jarray=60,deltay=1.0)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=30,imax=220,ioffset=-29,iarray=191,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 100.0E
c
cc      parameter (imin=30,imax=520,ioffset=-29,iarray=491,deltax=0.5)
cc      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)

c**   Parameter statement for 0.5 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 1.0W to 100.0E
c
      parameter (imin=2,imax=520,ioffset=-1,iarray=519,deltax=0.5)
      parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5)
c
c**   Parameter statement for 0.1 degree resolution for a total area of:
c**      latitude:   1.0N to  60.0N   longitude: 15.0W to 110.0W
c
cc      parameter (imin=150,imax=1100,ioffset=-149,iarray=951,deltax=0.1)
cc      parameter (jmin= 10,jmax= 600,joffset=  -9,jarray=591,deltay=0.1)
c
      integer numits, deltat, deltath
c
cc      parameter ( numits = 10, deltat = 12 )
      parameter ( numits = 20, deltat = 6 )
      parameter ( numitsh = 10, deltath = 12 )
c
c      common /final/ parray( iarray, jarray, 0:numits, 3 ),
c     &               parrayi( iarray, jarray, 0:numits, 3 )
     
      common /final/ parray( iarray, jarray, 0:numits, 3 ),
     &               parrayi( iarray, jarray, 0:numits, 3 ),
     &               parrayh( iarray, jarray, 0:numitsh, 3 )
     
c
      real       temparr(iarray,jarray)
      integer    lsec0(2),lsec1(13),igds(5),kpdt(200),kgdt(200)
      integer    idrt(200)
      integer    itype(3), itype2(3), end_plot, ngribunit, end_time
      logical(1) lb(iarray,jarray)
      character*1 cgrib(iarray*jarray*4)
c
      character*10 ymdh, ymdh_end
c
      data itype / 17491, 25722, 32924 /
      data itype2 / 198, 199, 200 /
      data ngribunit / 71 /
c
cc      lb = .TRUE.
      lcgrib = iarray*jarray*4
c
c**   Parse the initial dtg
c
      read (ymdh(1:4), '(i4)') iyyyy
      read (ymdh(1:2), '(i2)') icc
      read (ymdh(3:4), '(i2)') iyy
      read (ymdh(5:6), '(i2)') imm
      read (ymdh(7:8), '(i2)') idd
      read (ymdh(9:10),'(i2)') ihh
c
      icent = icc + 1
c
      print *,'iyyyy= ',iyyyy,'  imm= ',imm,' idd= ',idd,' ihh= ',ihh
c
      nxstart = (nint (360.0 - (float(iarray) * deltax) + 
     &                        (float(ioffset) * deltax))) * 1000000
      nyend   = (nint (float(jmax) * deltay)) * 1000000
c
      nxend   = nxstart + (nint (float(iarray-1) * deltax)) * 1000000
      nystart = nyend   - (nint (float(jarray-1) * deltay)) * 1000000
c
      ndeltax = nint (deltax * 1000000.0)
      ndeltay = nint (deltay * 1000000.0)
c
      kf = iarray*jarray
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 0
c
      lsec0(1) = 0              ! Discipline ( 0 = Meteorological )
      lsec0(2) = 2              ! GRIB Version
c
CSAG-------------------------------------------- 
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 1
c
      lsec1(1) = 7              ! Originating Center (TPC/NHC)
      lsec1(2) = 10             ! Sub-Center
      lsec1(3) = 2              ! Master Tables Version Number
      lsec1(4) = 1              ! Local Tables Version Number
      lsec1(5) = 0              ! Reference time is "Start of forecast"
cc      lsec1(5) = 1              ! Reference time is "Start of forecast"
      lsec1(6) = iyyyy          ! = initial yyyy of forecast
      lsec1(7) = imm            ! = initial   mm of forecast
      lsec1(8) = idd            ! = initial   dd of forecast
      lsec1(9) = ihh            ! = initial   hh of forecast
      lsec1(10) = 0             ! = initial minutes of forecast
      lsec1(11) = 0             ! = initial seconds of forecast
      lsec1(12) = 0             ! Status of data = Production
      lsec1(13) = 1             ! Type of data = Forecast
c
CSAG-------------------------------------------- 
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 3 - GDS
c
      igds(1) = 0               ! Grid definition in GDT
      igds(2) = kf              ! Number of grid points  
      igds(3) = 0               ! No additional point definitions
      igds(4) = 0               ! No optional List
      igds(5) = 0               ! Grid Definition Template No. (0 = lat/lon)  
c
      kgdt(1) = 6               ! Predefined Radius of earth
      kgdt(2) = 0               !
      kgdt(3) = 0               !
      kgdt(4) = 0               !
      kgdt(5) = 0               !
      kgdt(6) = 0               !
      kgdt(7) = 0               !
      kgdt(8) = iarray          ! = # pts in x-direction
      kgdt(9) = jarray          ! = # pts in y-direction
      kgdt(10) = 0              !
      kgdt(11) = 0              !
      kgdt(12) = nystart        ! = Latitude of first point (*1000000)
      kgdt(13) = nxstart        ! = Longitude of first point (*1000000)
      kgdt(14) = 48             ! Resolution flags
      kgdt(15) = nyend          ! = Latitude of last point (*1000000)
      kgdt(16) = nxend          ! = Longitude of last point (*1000000)
      kgdt(17) = ndeltax        ! = x-increment in degrees * 1000000
      kgdt(18) = ndeltay        ! = y-increment in degrees * 1000000
      kgdt(19) = 64             ! Scanning Mode
c
CSAG-------------------------------------------- 
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 4 - PDS 
c
c**   Defined later in the subroutine
c
cc      kpdtn = 9                 ! Product Definition Template Number
      kpdt=0
      kpdt(1) = 2               ! Parameter Category
      kpdt(2) = 1               ! Parameter Number
      kpdt(3) = 2               ! Forecast
      kpdt(4) = 0
      kpdt(5) = 30              ! Generating process
      kpdt(8) = 1               ! Time unit = hour
c
c**   If the incremental fields are ever defined kpdt(9) will be 
c**     defined for the beginning of each increment
c
cc      kpdt(9) = 0               ! Beginning hour for accumulation
      kpdt(10) = 103            ! 1st level = surface
      kpdt(12) = 10             ! height in meters above the ground 
      kpdt(13) = 255            ! 2nd level = missing
      kpdt(16) = 255            ! Probability number
      kpdt(17) = 255            ! total probabilities
      kpdt(18) = 1
      kpdt(19) = 0
      kpdt(20) = 0
      kpdt(21) = 3              ! scale factor for threshold
c
c**   These will be defined later in the subroutine
c**   When kpdtn = 5, only kpdt(1) through (22) is used
c
cc      kpdt(22) = (defined in code)  ! thresholds
cc      kpdt(23) = iyyyy_end      ! fcst end year
cc      kpdt(24) = imm_end        ! fcst end month
cc      kpdt(25) = idd_end        ! fcst end day
cc      kpdt(26) = ihh_end        ! fcst end hour
c
      kpdt(27) = 0              ! fcst end minutes
      kpdt(28) = 0              ! fcst end seconds
c
      kpdt(29) = 1              ! number of time ranges 
      kpdt(30) = 0              ! number missing
      kpdt(31) = 1              ! indicates accumulation
      kpdt(32) = 2              ! type of increment
      kpdt(33) = 1              ! time range unit (hours)
cc      kpdt(34) = (defined in code)  ! length of time range
      kpdt(35) = 1              ! time unit of increment
      kpdt(36) = 0              ! length of increment
c
CSAG-------------------------------------------- 
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 6
c
      ibmap=255
c
CSAG-------------------------------------------- 
c
CSAG-------------------------------------------- 
c
C**   Set info for GRIB2 Section 5 and 7 (reset in the loop, below)
c
cc      idrtn = 0                ! Data Representation Template Number
cc      idrt(1) = 0              ! Reference value - Set by GRIB2 encoder
cc      idrt(2) = 0              ! Binary scale factor
cc      idrt(3) = 1              ! Decimal scale factor
cc      idrt(4) = 0              ! Number of bits - Set by encoder, if 0
cc      idrt(5) = 0              ! Floating point data
c
CSAG-------------------------------------------- 
c
c**   Now run through for each forecast hour and call the NCEP
c**   w3lib routine "putgb", which packs the data into 1 output
c**   GRIB record and writes it to unit "ngribunit".
c
       
      call open_grib_files ( ngribunit, iret )
c
      if (iret /= 0) then
         print '(/,a50,2i4,/)', '!!! ERROR: opening unit, rc= ', 
     &        ngribunit, iret
         stop 99
      endif
c
cc      open  ( 32, file='grib2_input_array', status='unknown' )
c
      do n = 0, numits
c
         kpdtn = 9
c
cc         if ( n .eq. 0 ) kpdtn = 5
c
cc            kpdt(9) = n*deltat
c
         end_time = n*deltat
c
c**   Parse the initial dtg
c
         call dtgmod ( ymdh, end_time, ymdh_end, istat )
         if ( istat .ne. 0 ) stop 'Bad ymdh_end calculated!'
c
         read (ymdh_end(1: 4),'(i4)') iyyyy_end
         read (ymdh_end(5: 6),'(i2)') imm_end
         read (ymdh_end(7: 8),'(i2)') idd_end
         read (ymdh_end(9:10),'(i2)') ihh_end
c
         print *,'iyyyy_end= ',iyyyy_end,'  imm_end= ',imm_end,
     &        ' idd_end= ',idd_end,' ihh_end= ',ihh_end
c
         kpdt(23) = iyyyy_end   ! fcst end year
         kpdt(24) = imm_end     ! fcst end month
         kpdt(25) = idd_end     ! fcst end day
         kpdt(26) = ihh_end     ! fcst end hour
c
c**   Write out the 34-knot, 50-knot and 64-knot cummulative probabilities
c
         kpdt(9)  = 0
         kpdt(34) = n*deltat
         kpdt(32) = 192    ! type of increment - identifies cumulative probabilities (lauer/boyer)
c
         do m = 1, 3
c
            do j = 1, jarray
               do i = 1, iarray
c
c**   Flip the probability array in the x-direction
c
                  iflip = iarray - i + 1
                  temparr( iflip, j ) = parray( i, j, n, m )
c
c**   Set the bit map for valid data
c
                  lb( iflip, j ) = .true.
c  		
               enddo
            enddo
c
cc            write ( 32, '( 2i6 )' ) kpdt(34), m
cc            write ( 32, '(119f5.1)' ) 
cc     &           ((temparr( i, j ), j = 1, jarray), i = 1, iarray)
c
            kpdt(22) = itype(m)   ! wind speed threshold in m/s
c
CSAG-------------------------------------------- 
c
           call gribcreate (cgrib,lcgrib,lsec0,lsec1,iret)
            print *, 'gribcreate return code for grib2out: iret= ', iret
c
           call addgrid (cgrib,lcgrib,igds,kgdt,200,idummy,1,iret)
            print *, 'addgrid return code for grib2out: iret= ', iret
c
C**   Reset info for GRIB2 Section 5 and 7
c
            idrtn = 0           ! Data Representation Template Number
            idrt(1) = 0         ! Reference value - Set by GRIB2 encoder
            idrt(2) = 0         ! Binary scale factor
            idrt(3) = 1         ! Decimal scale factor
            idrt(4) = 0         ! Number of bits - Set by encoder, if 0
            idrt(5) = 0         ! Floating point data
c
            call addfield (cgrib,lcgrib,kpdtn,kpdt,200,idummy,0,idrtn,
     &              idrt,200,temparr,kf,ibmap,lb,iret)
            print *, 'addfield return code for grib2out: iret= ', iret
c
            call gribend (cgrib,lcgrib,lencgrib,iret)
            print *, 'gribend return code for grib2out: iret= ', iret
c
            call wryte (ngribunit,lencgrib,cgrib)
c
CSAG-------------------------------------------- 
c
         enddo
c
c**   Write out the 34-knot, 50-knot and 64-knot incremental probabilities
c
         if ( n .eq. 0 ) then
            kpdt(9)  = 0
            kpdt(34) = 0
         else
            kpdt(9)  = (n - 1)*deltat
            kpdt(34) = deltat
         endif
         kpdt(32) = 2 ! type of increment - identifies incremental (Lauer/Boyer)

c
         do m = 1, 3
c
            do j = 1, jarray
               do i = 1, iarray
c
c**   Flip the probability array in the x-direction
c
                  iflip = iarray - i + 1
                  temparr( iflip, j ) = parrayi( i, j, n, m )
c
c**   Set the bit map for valid data
c
                  lb( iflip, j ) = .true.
c  		
               enddo
            enddo
c
cc            write ( 32, '( 2i6 )' ) kpdt(34), m
cc            write ( 32, '(119f5.1)' ) 
cc     &           ((temparr( i, j ), j = 1, jarray), i = 1, iarray)
c
            kpdt(22) = itype(m)   ! wind speed threshold in m/s
c
CSAG-------------------------------------------- 
c
            call gribcreate (cgrib,lcgrib,lsec0,lsec1,iret)
            print *, 'gribcreate return code for grib2out: iret= ', iret
c
            call addgrid (cgrib,lcgrib,igds,kgdt,200,idummy,1,iret)
            print *, 'addgrid return code for grib2out: iret= ', iret
c
C**   Reset info for GRIB2 Section 5 and 7
c
            idrtn = 0           ! Data Representation Template Number
            idrt(1) = 0         ! Reference value - Set by GRIB2 encoder
            idrt(2) = 0         ! Binary scale factor
            idrt(3) = 1         ! Decimal scale factor
            idrt(4) = 0         ! Number of bits - Set by encoder, if 0
            idrt(5) = 0         ! Floating point data
c
            call addfield (cgrib,lcgrib,kpdtn,kpdt,200,idummy,0,idrtn,
     &              idrt,200,temparr,kf,ibmap,lb,iret)
            print *, 'addfield return code for grib2out: iret= ', iret
c
            call gribend (cgrib,lcgrib,lencgrib,iret)
            print *, 'gribend return code for grib2out: iret= ', iret
c
            call wryte (ngribunit,lencgrib,cgrib)
c
CSAG-------------------------------------------- 
c
         enddo
c
      enddo
      
      


c
c**   Write out the 34-knot, 50-knot and 64-knot incremental probabilities at 6 hourly increments, formatted for Pablo
c
      do n = 1, numits

         end_time = n*deltat
	 
	 call dtgmod ( ymdh, end_time, ymdh_end, istat )
         if ( istat .ne. 0 ) stop 'Bad ymdh_end calculated!'
c
         read (ymdh_end(1: 4),'(i4)') iyyyy_end
         read (ymdh_end(5: 6),'(i2)') imm_end
         read (ymdh_end(7: 8),'(i2)') idd_end
         read (ymdh_end(9:10),'(i2)') ihh_end
c
         print *,'iyyyy_end= ',iyyyy_end,'  imm_end= ',imm_end,
     &        ' idd_end= ',idd_end,' ihh_end= ',ihh_end
c
         kpdt(23) = iyyyy_end   ! fcst end year
         kpdt(24) = imm_end     ! fcst end month
         kpdt(25) = idd_end     ! fcst end day
         kpdt(26) = ihh_end     ! fcst end hour

        kpdt(9)  = (n - 1)*deltat
        kpdt(34) = deltat
        
         kpdt(32) = 2 ! type of increment - identifies incremental (Lauer/Boyer)

c
         do m = 1, 3
            do j = 1, jarray
               do i = 1, iarray
c
c**   Flip the probability array in the x-direction
c
                  iflip = iarray - i + 1
                  temparr( iflip, j ) = parrayi( i, j, n, m )
c
c**   Set the bit map for valid data
c
                  lb( iflip, j ) = .true.
c  		
               enddo
            enddo
c
cc            write ( 32, '( 2i6 )' ) kpdt(34), m
cc            write ( 32, '(119f5.1)' ) 
cc     &           ((temparr( i, j ), j = 1, jarray), i = 1, iarray)
c
            kpdt(22) = itype(m)   ! wind speed threshold in m/s
            kpdt(2) = itype2(m)      ! Parameter Numbers are 198,199,200 for incrementals
c
CSAG-------------------------------------------- 
c
            call gribcreate (cgrib,lcgrib,lsec0,lsec1,iret)
            print *, 'gribcreate return code for grib2out: iret= ', iret
c
            call addgrid (cgrib,lcgrib,igds,kgdt,200,idummy,1,iret)
            print *, 'addgrid return code for grib2out: iret= ', iret
c
C**   Reset info for GRIB2 Section 5 and 7
c
            idrtn = 0           ! Data Representation Template Number
            idrt(1) = 0         ! Reference value - Set by GRIB2 encoder
            idrt(2) = 0         ! Binary scale factor
            idrt(3) = 1         ! Decimal scale factor
            idrt(4) = 0         ! Number of bits - Set by encoder, if 0
            idrt(5) = 0         ! Floating point data
c
            call addfield (cgrib,lcgrib,kpdtn,kpdt,200,idummy,0,idrtn,
     &              idrt,200,temparr,kf,ibmap,lb,iret)
            print *, 'addfield return code for grib2out: iret= ', iret
c
            call gribend (cgrib,lcgrib,lencgrib,iret)
            print *, 'gribend return code for grib2out: iret= ', iret
c
            call wryte (ngribunit,lencgrib,cgrib)
c
CSAG-------------------------------------------- 
c
         enddo
c
      enddo
      
      
c
cc      close ( 32 )
c
      call baclose ( ngribunit, iret )
      print *, 'baclose return code for grib2out: iret= ', iret
c
      return
      end