      program degrib2model_driver
c
      implicit none
c

      integer nx,ny,nz,kpds(200),datsav(9),n         !Degrib grid dimensions
	integer filemax,noct,I
	parameter (filemax=8)
	integer octnum(filemax),jm1,jp1,ip1,ip2,ip3,ip4,latcov,loncov,
     +					readin
c
      character(LEN=255):: gribfile,outdir
c
      real esat,es
c
      common /estab/esat(15000:45000),es(15000:45000)
c_______________________________________________________________________________
c
c *** Initialize tables.
c
      call es_ini
c
c *** Read file names.
c
Cmp
Cmp	all in same file....read in desired octant numbers????
Cmp
      read(5,'(a)') gribfile
      read(5,'(a)') outdir

	do I=1,filemax+1
C	write(6,*) 'reading for I= ', I
	read(5,23) readin
	if (readin .eq. -9) then
	noct=I-1
	goto 1103
	endif
	octnum(I)=readin
	enddo

  23	format(I2)
 1103	continue
c

	write(6,*) 'want to work with the following octants: '
	do I=1,noct
	write(6,*) octnum(I)
	enddo

	
	n=index(gribfile,' ')-1

	IP1=0
	IP2=0
	IP3=0
	IP4=0

	JP1=0
	JM1=0

	do I=1,noct
C	CALL GET_GDS(gribfile(1:n),octnum(I),datsav,kpds)

	if (octnum(I) .le. 40) JP1=1
	if (octnum(I) .gt. 40) JM1=1
	if (octnum(I) .eq. 37 .or. octnum(I) .eq. 41) IP1=1
	if (octnum(I) .eq. 38 .or. octnum(I) .eq. 42) IP2=1
	if (octnum(I) .eq. 39 .or. octnum(I) .eq. 43) IP3=1
	if (octnum(I) .eq. 40 .or. octnum(I) .eq. 44) IP4=1

C	datsav(3) SW lat
C	datsav(4) SW lon

	enddo

	latcov=JM1+JP1
	loncov=(IP1+IP2+IP3+IP4)

	write(6,*) 'octants selected are ', loncov, ' in longitude '
	write(6,*) 'and ', latcov, ' in latitude'
C	write(6,*) 'calculated NX x NY is : ', loncov*73, latcov*73

        nz=11
C	nx=loncov*73
C	ny=latcov*73
	nx=loncov*73-(loncov-1)
	ny=latcov*73-(latcov-1)

	write(6,*) 'nx= ny= ', nx,ny
C	stop

      call degrib2model(gribfile,outdir,octnum,noct,nx,ny,nz)

c
      end
c
c===============================================================================
c
      subroutine degrib2model(gribfile,outdir,octnum,noct,nx,ny,nz)
c
      implicit none
c
      real    cp,kappa,dlat,dlon,sum,avg
      parameter(cp=1004.,kappa=287./1004.)
c
      integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld
     .         ,iyear,imonth,iday,ijulday,julday_laps,ifcsthr
     .	       ,count,jinput,noct,octnum(noct),swlatsv,swlonsv
      integer swlat(noct),swlon(noct),III,datsav(9),II
	integer londif,latdif,istart,jstart,ival,jval,kpds(200)

	integer KGDS(200)
c
      real ht(nx,ny,nz),htin(73,73,nz)   !Isobaric heights (m)
     .      ,tp(nx,ny,nz),tpin(73,73,nz)   !Isobaric temps (K)
     .      ,th(nx,ny,nz),thin(73,73,nz)   !Isobaric theta (K)
     .      ,uw(nx,ny,nz),uwin(73,73,nz)   !Isobaric u-wind (m/s)
     .      ,vw(nx,ny,nz),vwin(73,73,nz)   !Isobaric v-wind (m/s)
     .      ,rh(nx,ny,nz),rhin(73,73,nz)   !Isobaric rh,mr (%,kg/kg)
     .      ,ex(nx,ny,nz)             !Isobaric Exner
     .      ,pr(nz),pri(nz)            !Isobaric pressures (mb)
     .      ,lat1,lat2,lon0,sw(2),ne(2)
     .      ,xe,mrsat,esat,es
     .	    ,avnlevs(11)
     .      ,temp(nx,ny)
     .	    ,slp(nx,ny,12)
c
      character(LEN=255):: gribfile,outdir,outfile,gdsfile
      character(LEN=2)::   gproj
      character(LEN=11)::  atime
      character(LEN=8)::   model

       data AVNLEVS/1000,850,700,500,400,300,250,200,150,100,70/
c
      common /estab/esat(15000:45000),es(15000:45000)
c___________________________________________________________________________

Cmp Loop over whole subroutine
	n=index(gribfile,' ')-1
        do III=1,noct
	CALL GET_GDS(gribfile(1:n),octnum(III),datsav,kpds)
	swlat(III)=datsav(3)
	swlon(III)=datsav(4)
	enddo

        do I=1,noct
        swlatsv=99999
        if (swlat(I) .eq. -90000) then
        swlatsv=-90000
        goto 97
        endif
        enddo

C       lat not 90S, so must be 0
        swlatsv=0

  97    continue


Cmp	logic here is BAD!  
C
C	Need to be able to wrap from Grid 40 (120W at west) to grid
C	37 (30W at west).  
C
        do I=1,noct
        if (swlon(I) .eq. -30000) then
        swlonsv=-30000
        goto 98
        endif
        enddo

        do I=1,noct
        if (swlon(I) .eq. 60000) then
        swlonsv=60000
        goto 98
        endif
        enddo

        do I=1,noct
        if (swlon(I) .eq. 150000) then
        swlonsv=150000
        goto 98
        endif
        enddo

        swlonsv=-120000

  98    continue

C
C
C	at this point have initial guess of "westernmost" point, 
C	but not sure if it is the "left" most point on the final grid

	if (swlonsv .eq. -30000) then
  	  do III=1,noct

	    if (swlon(III) .eq. -120000) then
	      write(6,*) 'special case...to 120W'
	      swlonsv=-120000
	      do II=1,noct
	        if (swlon(II) .eq. 150000) then
	          write(6,*) 'special case...to 150E'
	          swlonsv=150000
		endif
	      enddo
            endif

          enddo
	endif

	do III=1,noct
	jinput=octnum(III)
c
c *** Fill pressure levels.
c
C	write(6,*) 'using hardwired stuff! ', nx,ny,nz
      do k=1,nz
	pr(k)=AVNLEVS(k)
C		write(6,*) 'pressure value of ', pr(k)
	      enddo
c 
c *** Read in degrib data.
c
Cmp
Cmp	"thinned" WAFS grids have 73 rows, and when filled are 73 X 73
Cmp	However, in thinned status each octant has 3447 gridpoints
Cmp
      n=index(gribfile,' ')-1
      call read_degrib(gribfile(1:n)
     .     ,3447,nz,pr,htin,tpin,rhin,uwin,vwin,jinput,atime)

C
C Find destination 

	londif=swlon(III)-swlonsv
	latdif=swlat(III)-swlatsv
	if (londif .lt. 0) londif=londif+360000
	if ( londif .eq. 0) ISTART=1
C	if ( londif.eq.90000) ISTART=74
	if ( londif.eq.90000) ISTART=73
C	if ( londif.eq.180000) ISTART=147
	if ( londif.eq.180000) ISTART=145
C	if ( londif.eq.270000) ISTART=146+73
	if ( londif.eq.270000) ISTART=145+72
	if ( latdif .eq. 0) JSTART=1
C	if ( latdif .eq. 90000) JSTART=74
	if ( latdif .eq. 90000) JSTART=73
C
C
C Put the "in" arrays into their final places for the isobaric fields
C

	do  K=1,NZ
	do  J=JSTART,JSTART+73-1
	do  I=ISTART,ISTART+73-1

	ival=(I-istart)+1
        jval=(J-jstart)+1

        ht(I,J,K)=htin(ival,jval,K)
        tp(I,J,K)=tpin(ival,jval,K)
        rh(I,J,K)=rhin(ival,jval,K)
        uw(I,J,K)=uwin(ival,jval,K)
        vw(I,J,K)=vwin(ival,jval,K)
	
	enddo
	enddo
	enddo


c
c *** Check for any missing data.
c

	if (III .eq. NOCT) THEN
      do k=1,nz
	  if (ht(1,1,k) .eq. -99999.) then
	     print *,'Height data missing at level: ',pr(k)
	     stop
	  elseif (tp(1,1,k) .eq. -99999.) then
	     print *,'Temperature data missing at level: ',pr(k)
	     stop
	  elseif (rh(1,1,k) .eq. -99999.) then
	     print *,'RH data missing at level: ',pr(k)
	     print *,'Calling RH patch.'
	     call rh_fix(nx,ny,nz,rh,pr)
	  elseif (uw(1,1,k) .eq. -99999.) then
	     print *,'U-wind data missing at level: ',pr(k)
	     stop
	  elseif (vw(1,1,k) .eq. -99999.) then
	     print *,'V-wind data missing at level: ',pr(k)
	     stop
	  endif
      enddo

c
c *** Convert 3d temp to theta.
c *** Compute Exner function.
c *** Convert 3d rh to mr.
c
      do k=1,nz
	 pri(k)=1./pr(k)
      enddo
c
      do k=1,nz
      do j=1,ny
      do i=1,nx
	 th(i,j,k)=tp(i,j,k)*(1000.*pri(k))**kappa
	 ex(i,j,k)=cp*tp(i,j,k)/th(i,j,k)
	 it=tp(i,j,k)*100
	 it=min(45000,max(15000,it))
	 xe=esat(it)
	 mrsat=0.00622*xe/(pr(k)-xe)
	 rh(i,j,k)=rh(i,j,k)*mrsat
      enddo
      enddo
      enddo


Cmp
Cmp for some reason AVN data appears to be "upside down" in the N-S
Cmp sense.  Flip the arrays of data.
Cmp

C      do k=1,nz

C      do j=1,ny
C      do i=1,nx
C	temp(i,j)=th(i,j,k)
C      enddo
C      enddo

C      do J=1,ny
C      do i=1,nx
C	th(i,j,k)=temp(i,ny-J+1)
C      enddo
C      enddo

C======================================

C	      do j=1,ny
C	      do i=1,nx
C		temp(i,j)=uw(i,j,k)
C	      enddo
C	      enddo

C	      do J=1,ny
C	      do i=1,nx
C		uw(i,j,k)=temp(i,ny-J+1)
C	      enddo
C	      enddo

C======================================

C	      do j=1,ny
C	      do i=1,nx
C		temp(i,j)=vw(i,j,k)
C	      enddo
C	      enddo

C	      do J=1,ny
C	      do i=1,nx
C		vw(i,j,k)=temp(i,ny-J+1)
C	      enddo
C	      enddo

C======================================

C	      do j=1,ny
C	      do i=1,nx
C		temp(i,j)=ht(i,j,k)
C	      enddo
C	      enddo

C		if (k.eq.1) write(6,*) 'flipping isobaric data - '
C	      do J=1,ny
C	      do i=1,nx
C		ht(i,j,k)=temp(i,ny-J+1)
C	      enddo
C	      enddo

C======================================

C	      do j=1,ny
C	      do i=1,nx
C		temp(i,j)=rh(i,j,k)
C	      enddo
C	      enddo

C	      do J=1,ny
C	      do i=1,nx
C		rh(i,j,k)=temp(i,ny-J+1)
C	      enddo
C	      enddo

C======================================

C	      do j=1,ny
C	      do i=1,nx
C		temp(i,j)=ex(i,j,k)
C	      enddo
C	      enddo

C	      do J=1,ny
C	      do i=1,nx
C		ex(i,j,k)=temp(i,ny-J+1)
C	      enddo
C	      enddo
C-----------------------------------------
C      enddo

c
C      print *,'ua :',ht(nx/2,ny/2,nz),th(nx/2,ny/2,nz),rh(nx/2,ny/2,nz)
C     .       ,uw(nx/2,ny/2,nz),vw(nx/2,ny/2,nz)
C      print *,'ua :',ht(nx/2,ny/2,1),th(nx/2,ny/2,1),rh(nx/2,ny/2,1)
C    .       ,uw(nx/2,ny/2,1),vw(nx/2,ny/2,1)
c
c *** Create output file name.
c
      n=index(outdir,' ')-1

      model='ETA_wafs'
      outfile=outdir(1:n)//'/'//atime//'.'//model
      n=index(outfile,' ')-1
      print *,model,' data ---> ',outfile(1:n)
      open(1,file=outfile(1:n),status='unknown',form='unformatted')
c
c *** Write header stuff.

c
      ip=1
      jp=1
Cmp      nsfcfld will be 12 although 4 of these are set to missing for AVN...
      nsfcfld=12
Cmp	is the AVN 'LL' or 'CE'??????
      gproj='LL'
!      write(1) nx,ny,nz,nx,ny,ip,jp,nsfcfld,gproj
	write(1) nz

C new stuff for GDS file
C
        n=index(gribfile,' ')-1
        write(6,*) 'calling get_fullgds with ',
     +    gribfile(1:n)
        call get_fullgds(gribfile(1:n),nx,ny,kgds)
        write(6,*) 'back from get_fullgds'


Cmoved        n=index(outdir,' ')-1
C        gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model
C
C        n=index(gdsfile,' ')-1
C        write(6,*) 'GDS written to ', gdsfile(1:n)
C        open(unit=14,file=gdsfile(1:n),form='unformatted',
C     +  access='sequential')
C
C        write(6,*) 'writing kgds ', (kgds(I),i=1,14)
C        write(14) kgds
C
Cmoved        close(14)
C
C end new stuff

C
C
C   this header info is octant-dependent...also needs work for patched
C	together octants
C
C      lat1=-90.0
C      lat2=0.0
C      lon0=0.0
C      sw(1)=-90
C      sw(2)=0.
C      ne(1)=90.
C      ne(2)=359.
	lat1=swlatsv/1000.
	lat2=0.
	lon0=swlonsv/1000.
	sw(1)=swlatsv/1000.
	sw(2)=swlonsv/1000.

Cmp	We know the SW corner coordinates, we know that each octant is 90 deg
Cmp	by 90 deg, and we know how many octants we are dealing with.  Calc
Cmp	the NE corner

	write(6,*) 'calculating the corner points assuming 
     +			this is SW corner: '
	write(6,*) swlatsv,swlonsv
	write(6,*) 'and nx and ny is this: ', nx,ny

	ne(1)= ( swlatsv + int(ny/73.+0.5)*90000 )/1000.
	ne(2)= ( swlonsv + int(nx/73.+0.5)*90000 )/1000.

	if (ne(2) .gt.  180.) ne(2)=ne(2)-360.
	
	write(6,*) 'using these values for the corners!!!!'
	write(6,*) 'SW ', sw(1),sw(2)
	write(6,*) 'NE ', ne(1),ne(2)
	dlat=1.25
	dlon=1.25
Cmp      write(1) nx,ny,nz,lat1,lat2,lon0,sw,ne
!	write(1) nx,ny,nz,lat1,lon0,dlat,dlon
c
c *** Write isobaric upper air data.
c

Cmp	write a dummy value for WAFS slp field to prevent
Cmp	read errors.

	slp=-9999.

      write(1) uw
      write(1) vw
C      write(1) th
      write(1) ht
      write(1) rh
C      write(1) ex
	write(1) pr
      write(1) slp
Cmp	

        n=index(outdir,' ')-1
        gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model

        n=index(gdsfile,' ')-1
        write(6,*) 'GDS written to ', gdsfile(1:n)
        open(unit=14,file=gdsfile(1:n),form='unformatted',
     +  access='sequential')

	kgds(2)=nx
	kgds(3)=ny
	kgds(4)=sw(1)*1000
	kgds(5)=sw(2)*1000
	kgds(7)=ne(1)*1000
	kgds(8)=ne(2)*1000

        write(6,*) 'writing kgds ', (kgds(I),i=1,14)
        write(14) kgds

        close(14)
C	write(6,*) 'writing these values to the ETA_wafs file...'
	do k=1,nz
	write(6,*) 'level, value ', k, ht(nx/2,ny/2,k)
C	write(6,*) 'level, value ', k, th(nx/2,ny/2,k)
C	write(6,*) 'level, value ', k, uw(nx/2,ny/2,k)
C	write(6,*) 'level, value ', k, vw(nx/2,ny/2,k)
C	write(6,*) 'level, value ', k, ex(nx/2,ny/2,k)
	enddo
      close(1)
c
	endif
	enddo
C	end loop III

      return
      end
c
c===============================================================================
c
      subroutine read_degrib(etafile,nxny,nz,pr
     .                      ,ht,tp,rh,uw,vw,jinput,atime)
c
      implicit none
c
      integer nxny,nz,i,nx,ny
	parameter(nx=73,ny=73)
c
      real pr(nz),ht(nx,ny,nz),tp(nx,ny,nz)
     .      ,rh(nx,ny,nz),uw(nx,ny,nz),vw(nx,ny,nz)
     .      ,dummy,waftmp(nxny)
c
      integer kgds(200),kpds(200),len,kerr
     .         ,lenpds,lenkgds,nwords,kpdsl
     .         ,ijulday,julday_laps,j,k,datsav(9)
     .	       ,IRETO,IRET1,JPDS(200),JGDS(200),KNUM,jinput
     .	       ,HEMFLAG

	logical BITMAP(nxny)
c
      character*(*) etafile
      character(LEN=1)::   pds(200)
      character(LEN=11)::   atime
c_______________________________________________________________________________
c
      len=index(etafile//' ',' ')-1
c
c
c *** Fill time stamp (assume all records are for same time).
c
	call get_gds(etafile(1:len),jinput,datsav,kpds)

C      ijulday=julday_laps(kpds(10),kpds(9),kpds(8))
	if (kpds(8) .eq. 100) kpds(8)=0
      write(atime,'(4(i2.2),i3.3)') 
     .   kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)

	write(6,*) 'atime= ', atime
c
c *** Fill a missing value flag into first space of each variable.
c
      do k=1,nz
	 ht(1,1,k)=-99999.
 	 tp(1,1,k)=-99999.
	 rh(1,1,k)=-99999.
	 uw(1,1,k)=-99999.
	 vw(1,1,k)=-99999.
      enddo

c
c *** Now put the data into the corresponding arrays.
c

	call baopen(11,etafile,IRETO)
	if (IRETO .ne. 0) write(6,*) 'BAOPEN TROUBLE!!!! ', IRETO

	jpds=-1
	jgds=-1

	jpds(3)=jinput
	
C 	ISOBARIC DATA
C
C
C **************************************************
C
C	NOTE::::NEED TO DISTINGUISH BETWEEN NH AND SH OCTANTS
C	AS THEY GIVE A DIFFERENT ARGUMENT TO W3FT33!!!!!!!!!!!!!!
C
C **************************************************
C
C

        jpds(6)=100

	if (jpds(3).le.40) HEMFLAG=1
	if (jpds(3).ge.41) HEMFLAG=-1

        do k=1,nz
        jpds(7)=nint(pr(k))

        jpds(5)=7
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,waftmp,IRET1)
	if (IRET1 .eq. 0) call w3ft33(waftmp,ht(1,1,k),HEMFLAG)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)

        jpds(5)=11
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,waftmp,IRET1)
	if (IRET1 .eq. 0) call w3ft33(waftmp,tp(1,1,k),HEMFLAG)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)

        jpds(5)=52
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,waftmp,IRET1)
	if (IRET1 .eq. 0) call w3ft33(waftmp,rh(1,1,k),HEMFLAG)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)

        jpds(5)=33
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,waftmp,IRET1)
	if (IRET1 .eq. 0) call w3ft33(waftmp,uw(1,1,k),HEMFLAG)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)

        jpds(5)=34
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,waftmp,IRET1)
	if (IRET1 .eq. 0) call w3ft33(waftmp,vw(1,1,k),HEMFLAG)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)


C        write(6,*) 'Z,T,Q,U,V ', ht(nx/2,ny/2,k),tp(nx/2,ny/2,k),
C     +  rh(nx/2,ny/2,k),uw(nx/2,ny/2,k),vw(nx/2,ny/2,k)

       enddo

c
c *** Normal finish.
c
1000  continue
      return
c
c *** Premature end of file.
c
1100  continue
      print *,'Premature end of file.'
      print *,'Abort...'
      stop
c
      end
c
c===============================================================================
c
      subroutine rh_fix(nx,ny,nz,rh,pr)
c
      implicit none
c
      integer nx,ny,nz,i,j,k,kk
c
      real rh(nx,ny,nz),pr(nz)
c_______________________________________________________________________________
c
c *** Fix bottom levels if necessary.
c
      if (rh(1,1,1) .eq. -99999.) then
	 do k=2,nz
	    if (rh(1,1,k) .ne. -99999.) then
	       do kk=k-1,1,-1
	       print *,'Copying',nint(pr(kk+1)),' mb to'
     .                , nint(pr(kk)),' mb.'
	       do j=1,ny
	       do i=1,nx
		  rh(i,j,kk)=rh(i,j,kk+1)
               enddo
               enddo
               enddo
	       goto 10
            endif
         enddo
	 print *,'RH patch did not work.'
	 stop
      endif
c
c *** Fix upper levels if necessary.
c
10    continue
      if (rh(1,1,nz) .eq. -99999.) then
	 do k=nz-1,1,-1
	    if (rh(1,1,k) .ne. -99999.) then
	       do kk=k+1,nz
	       print *,'Copying',nint(pr(kk-1)),' mb to'
     .                , nint(pr(kk)),' mb.'
	       do j=1,ny
	       do i=1,nx
		  rh(i,j,kk)=rh(i,j,kk-1)
               enddo
               enddo
               enddo
	       goto 20
            endif
         enddo      
      endif
c
20    continue
      do k=1,nz
	 if (rh(1,1,k) .eq. -99999.) then
	    print *,'RH patch did not work.'
	    stop
	 endif
      enddo
c
      return
      end
c
c===============================================================================
c
      subroutine es_ini
c
      common /estab/esat(15000:45000),es(15000:45000)
c
c *** Create tables of the saturation vapour pressure with up to
c        two decimal figures of accuraccy:
c
      do it=15000,45000
         t=it*0.01
         p1 = 11.344-0.0303998*t
         p2 = 3.49149-1302.8844/t
         c1 = 23.832241-5.02808*alog10(t)
         esat(it) = 10.**(c1-1.3816E-7*10.**p1+
     .               8.1328E-3*10.**p2-2949.076/t)
         es(it) = 610.78*exp(17.269*(t-273.16)/(t-35.86))
      enddo
c
      return
      end
c
c===============================================================================
c
      function julday_laps(day,month,year)
c
      implicit	none
c      
      integer	julday_laps
     .         ,day,month,year
     .         ,ndays(12),i
c
      data ndays/31,28,31,30,31,30,31,31,30,31,30,31/
c_______________________________________________________________________________
c
      julday_laps=0
      do i=1,month-1
         julday_laps=julday_laps+ndays(i)
      enddo
      julday_laps=julday_laps+day
      if (mod(year,4) .eq. 0 .and. month .gt. 2) 
     .   julday_laps=julday_laps+1
      return
c      
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       subroutine get_gds(etafile,jinput,gdsinfo,kpds)

        character*(*) etafile
        character(LEN=1):: pds(50)

Ctst        integer*4 kgds(200),kpds(50),len,kerr
        integer kgds(200),kpds(200),len,kerr
     .         ,lenpds,lenkgds,nwords,kpdsl
     .         ,ijulday,julday_laps,j,k,gdsinfo(9),jinput
     .         ,gdsav,IRETO,JGDS(200),JPDS(200)
	real tmp(6000)

	LOGICAL BITMAP(6000)

        nxny=73*73

        JPDS=-1
        JGDS=-1

        len=index(etafile//' ',' ')-1

        call baopen(11,etafile(1:len),IRETO)
C        write(6,*) 'BAOPEN in get_gds: ', IRETO

        if (IRETO .ne. 0) then
         print *,'Error opening unit=11, file name = ',etafile(1:len)
     .          ,' iostat = ',kerr
         stop
        endif

        jpds(5)=7
        jpds(6)=100
        jpds(7)=500
Cmp set octant
	jpds(3)=jinput
Cmp
C	write(6,*) 'calling getgb'
        call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,tmp,IRET1)
C	write(6,*) 'return from getgb ', IRET1


       gdsinfo(1)=KGDS(2)
       gdsinfo(2)=KGDS(3)
       gdsinfo(3)=KGDS(4)
       gdsinfo(4)=KGDS(5)
       gdsinfo(5)=KGDS(7)
       gdsinfo(6)=KGDS(8)
       gdsinfo(7)=KGDS(9)
       gdsinfo(8)=KGDS(12)
       gdsinfo(9)=KGDS(13)
Cmp       write(6,*) gdsinfo

        return
        print *,'GETGDS PROBLEM'
        stop

        end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE W3FT33(AIN,OUT,NSFLAG)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    W3FT33      THICKEN THINNED WAFS GRIB GRID 37-44
C   PRGMMR: RALPH PETTERSON  ORG: W/NMCXX    DATE: 94-11-13
C
C ABSTRACT: SUBROUTINE THICKENS ONE THINNED WAFS GRIB GRID TO A
C   REAL ARRAY OF 5329 NUMBERS (73,73) 1.25 DEGREE GRID.
C
C PROGRAM HISTORY LOG:
C   94-??-??  RALPH PETERSON
C   94-11-07  R.E.JONES        ADD DOC BLOCK, CHANGE CALL TO 3
C                              PARAMETERS. REPLACE COS WITH TABLE
C                              LOOKUP.
C   95-06-02  RALPH PETERSON   CHANGES TO CORRECT MISS-POSITION
C                              BETWEEN + OR - 8.75 N/S.
C   95-06-03  R.E.JONES        CHANGES SO 8 ROWS WITH 73 VALUES
C                              ARE NOT THICKENED, 10% FASTER.
C
C USAGE:    CALL W3FT33(AIN, OUT, NSFLAG)
C   INPUT ARGUMENT LIST:
C     AIN      - REAL 3447 WORD ARRAY WITH UNPACKED THINNED WAFS
C                GRIB TYPE 37-44.
C     NSFLAG   - INTEGER =  1  AIN IS WAFS GRIB GRID 37-40  N. HEMI.
C                        = -1  AIN IS WAFS GRIB GRID 41-44  S. HEMI.
C
C   OUTPUT ARGUMENT LIST:
C     OUT      - REAL (73,73) WORD ARRAY WITH THICKENED WAFS GRIB
C                GRID 37-44.
C
C REMARKS: THE POLE POINT FOR U AND V WIND COMPONENTS WILL HAVE ONLY
C   ONE POINT. IF YOU NEED THE POLE ROW CORRECTED SEE PAGE 9 SECTION
C   1 IN OFFICE NOTE 388. YOU NEED BOTH U AND V TO MAKE THE
C   CORRECTION.
C
C ATTRIBUTES:
C   LANGUAGE: SiliconGraphics 5.2 FORTRAN 77
C   MACHINE:  SiliconGraphics IRIS-4D/25, 35, INDIGO, Indy
C
C$$$
C
         PARAMETER (NX=73,NY=73)
         PARAMETER (NIN=3447)
C
         REAL      AIN(*)
         REAL      OUT(NX,NY)
C
         INTEGER   IPOINT(NX)
C
         SAVE
C
      DATA  IPOINT/
     & 73, 73, 73, 73, 73, 73, 73, 73, 72, 72, 72, 71, 71, 71, 70,
     & 70, 69, 69, 68, 67, 67, 66, 65, 65, 64, 63, 62, 61, 60, 60,
     & 59, 58, 57, 56, 55, 54, 52, 51, 50, 49, 48, 47, 45, 44, 43,
     & 42, 40, 39, 38, 36, 35, 33, 32, 30, 29, 28, 26, 25, 23, 22,
     & 20, 19, 17, 16, 14, 12, 11,  9,  8,  6,  5,  3,  2/
C
         NXM   = NX - 1
         FNXM  = FLOAT(NXM)
C
C        TEST FOR GRIDS (37-40)
C
         IF (NSFLAG.GT.0) THEN
C
C          DO NOT THICKEN 8 ROWS WITH 73 VALUES, MOVE DATA
C          TO OUT ARRAY. GRIDS (37-40) N.
C
           IS = 0
           DO J = 1,8
             DO I = 1,NX
               IS       = IS + 1
               OUT(I,J) = AIN(IS)
             END DO
           END DO
C
           IE = NX * 8
           DO J = 9,NY
             NPOINT   = IPOINT(J)
             IS       = IE + 1
             IE       = IS + NPOINT - 1
             DPTS     = (FLOAT(NPOINT)-1.) / FNXM
             PW       = 1.0
             PE       = PW + DPTS
             OUT(1,J) = AIN(IS)
             VALW     = AIN(IS)
             VALE     = AIN(IS+1)
             DVAL     = (VALE-VALW)
             DO I = 2,NXM
               WGHT     = PE -FLOAT(IFIX(PE))
               OUT(I,J) = VALW + WGHT * DVAL
               PW       = PE
               PE       = PE + DPTS
               IF (IFIX(PW).NE.IFIX(PE)) THEN
                 IS   = IS + 1
                 VALW = VALE
                 VALE = AIN(IS+1)
                 DVAL = (VALE - VALW)
               END IF
             END DO
             OUT(NX,J) = AIN(IE)
           END DO
C
         ELSE
C
C         DO NOT THICKEN 8 ROWS WITH 73 VALUES, MOVE DATA
C         TO OUT ARRAY. GRIDS (41-44) S.
C
           IS = NIN - (8 * NX)
           DO J = 66,NY
             DO I = 1,NX
               IS       = IS + 1
               OUT(I,J) = AIN(IS)
             END DO
           END DO
C
           IE = 0
           DO J = 1,65
             NPOINT   = IPOINT(74-J)
             IS       = IE + 1
             IE       = IS + NPOINT - 1
             DPTS     = (FLOAT(NPOINT)-1.) / FNXM
             PW       = 1.0
             PE       = PW + DPTS
             OUT(1,J) = AIN(IS)
             VALW     = AIN(IS)
             VALE     = AIN(IS+1)
             DVAL     = (VALE-VALW)
             DO I = 2,NXM
               WGHT     = PE -FLOAT(IFIX(PE))
               OUT(I,J) = VALW + WGHT * DVAL
               PW       = PE
               PE       = PE + DPTS
               IF (IFIX(PW).NE.IFIX(PE)) THEN
                 IS   = IS + 1
                 VALW = VALE
                 VALE = AIN(IS+1)
                 DVAL = (VALE - VALW)
               END IF
             END DO
             OUT(NX,J) = AIN(IE)
           END DO
         END IF
C
         RETURN
         END


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

                subroutine get_fullgds(etafile,nx,ny,kgds)

        character*(*) etafile
        character(LEN=1) pds(50)

        integer kgds(200),kpds(200),len,kerr,jpds(200)
     .         ,lenpds,lenkgds,nwords,kpdsl,jgds(200)
     .         ,j,k,KNUM,nx,ny

        logical bitmap(nx*ny)
        real tmp(nx*ny)

        write(6,*) 'inside get_fullgds...', etafile


        nxny=nx*ny

        len=index(etafile//' ',' ')-1

        jpds=-1
        jgds=-1

        jpds(5)=11
        jpds(6)=100
        jpds(7)=500

        write(6,*) 'calling getgb '
        write(6,*) 'jpds =  ', jpds
        write(6,*) 'jgds =  ', jgds

        call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,tmp,IRET1)
        write(6,*) 'return from getgb ', IRET1

        if (IRET1 .ne. 0) then
         print *,'Error  getting GDS in get_fullgds ', IRET1
         stop
        endif

        write(6,*) 'leaving get_fullgds'
        return
        end

