      program degrib2model_driver
c
      implicit none
c
      integer nx,ny,nz,kpds(200),datsav(9),n         !Degrib grid dimensions
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
      read(5,'(a)') gribfile
      read(5,'(a)') outdir
c
	n=index(gribfile,' ')-1
	write(6,*) 'first call of GET_GDS'
	CALL GET_GDS(gribfile(1:n),datsav,kpds)
      nx=datsav(1)
      ny=datsav(2)

	if (nx .eq.  360 .and. ny .eq. 181) then
        nz=26
	elseif (nx .eq. 144 .and. ny .eq. 73) then
	nz=17
	else
	write(6,*) 'what grid is this??? ', nx,ny
	write(6,*) ' '
	write(6,*) 'I do not know the vertical levels for this '
	write(6,*) 'grid...please add some code to dgeta2model.f'
	write(6,*) ' '
	STOP 999
	endif
c
Cold      call degrib2model(degrib_dir,degrib_file,outdir,nx,ny,nz)
      call degrib2model(gribfile,outdir,nx,ny,nz)

c
      end
c
c===============================================================================
c
Cold      subroutine degrib2model(degrib_dir,degrib_file,outdir,nx,ny,nz)
      subroutine degrib2model(gribfile,outdir,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,ifcsthr
     .	       ,count,kgds(200)
c
      real ht(nx,ny,nz)              !Isobaric heights (m)
     .      ,tp(nx,ny,nz)              !Isobaric temps (K)
     .      ,th(nx,ny,nz)              !Isobaric theta (K)
     .      ,uw(nx,ny,nz)              !Isobaric u-wind (m/s)
     .      ,vw(nx,ny,nz)              !Isobaric v-wind (m/s)
     .      ,rh(nx,ny,nz)              !Isobaric rh,mr (%,kg/kg)
     .      ,pr(nz),pri(nz)            !Isobaric pressures (mb)
     .      ,lat1,lat2,lon0,sw(2),ne(2)
     .      ,xe,mrsat,esat,es
     .	    ,avnlevs(26),reanlevs(17)
     .	    ,slp(nx,ny,12) 	  !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL
     .      ,temp(nx,ny)
				  !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC
				  !slp(i,j,5 & 6) are 0-10 cm STC and SMC
				  !slp(i,j,7 & 8) are 10-200 cm STC and SMC

	REAL, ALLOCATABLE:: htout(:,:,:),
     +			    uwout(:,:,:),vwout(:,:,:),rhout(:,:,:),
     +			   slpout(:,:,:)
c
      character*255 gribfile,outdir,outfile,gdsfile
      character*2   gproj
      character*12   atime
      character*7   model

       data AVNLEVS/1000,975,950,925,900,850,800,750,700,650,
     +		     600,550,500,450,400,350,300,250,200,150,
     +		     100, 70, 50, 30, 20, 10/

	data REANLEVS/1000,925,850,700,600,500,
     +		       400,300,250,200,150,100,
     +		        70, 50, 30, 20, 10/
c
      common /estab/esat(15000:45000),es(15000:45000)
c_______________________________________________________________________________
c
c *** Fill pressure levels.
c
Cmp      nsfcfld will be 12 although 4 of these are set to missing for AVN...
      nsfcfld=12
Cmp
	write(6,*) 'using hardwired stuff! ', nx,ny,nz
      do k=1,nz
	if (nz .eq. 26) then
  	  pr(k)=AVNLEVS(k)
	elseif (nz .eq. 17) then
	  pr(k)=reanlevs(k)
	endif
C	write(6,*) 'pr(K): ', pr(K)
      enddo
c 
c *** Read in degrib data.
c
      n=index(gribfile,' ')-1
      write(6,*) 'call read_degrib'
      call read_degrib(gribfile(1:n)
     .                ,nx*ny,nz,pr,ht,tp,rh,uw,vw,slp,atime)
	write(6,*) 'return read_degrib'
c
c *** Check for any missing data.
c
      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

Cmp
Cmp     Handle missing surface data....
Cmp

C
C	For reanalysis runs these data will not be used.  Just bogus
C	to hold a place.
C
	write(6,*) 'going to use nsfcfld= ', nsfcfld

C	If rean put bogus data into slp

	if (nz .eq. 17) then
	  do K=5,8
	do j=1,ny
	do i=1,nx
	    if (mod(k,2) .eq. 0) then
	      slp(i,j,k)=0.14
	    else
	      slp(i,j,k)=275.
	    endif
	enddo
	enddo
	  enddo
	endif

	do k=5,nsfcfld

	if (slp(1,1,k) .eq. -99999.) then
 	write(6,*) 'SURFACE DATA MISSING... ', K       
	if (k.eq.9.or.k.eq.11) then
	 write(6,*) 'filling soil temp data...'
	 do j=1,ny
	 do i=1,nx
	slp(i,j,k)=slp(i,j,7)
	 enddo
	 enddo

      elseif(k.eq.10.or.k.eq.12.) then

	write(6,*) 'filling soil moisture data...'

	do j=1,ny
	do i=1,nx
	 slp(i,j,k)=slp(i,j,8)
	enddo
	enddo
      endif
      endif
Cmp
	if (mod(k,2) .eq. 0) then
	do j=1,ny
	do i=1,nx
	if (slp(i,j,k) .eq. 0) then
	slp(i,j,k)=0.14
	endif
	enddo
	enddo

	else

	do j=1,ny
        do i=1,nx
        if (slp(i,j,k) .eq. 0) then
        slp(i,j,k)=273.15
        endif
        enddo
        enddo

	endif
Cmp
	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
C	 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.  If reanalysis data, skip this part
Cmp
	if (nz .eq. 17) goto 969
      do k=1,nz

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Cmp	generalize eventually!
	do k=1,nsfcfld

C	write(6,*) 'flipping surface fields... ' 
	
      do j=1,ny
      do i=1,nx
	temp(i,j)=slp(i,j,k)
      enddo
      enddo

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

	enddo

  969	continue

CC
CC
CC 	Add bogus column if reanalysis data
CC
CC

	if (nx .eq. 144) then
	ALLOCATE(HTOUT(145,ny,nz),UWOUT(145,ny,nz),VWOUT(145,ny,nz))
	ALLOCATE(RHOUT(145,ny,nz),SLPOUT(145,ny,nsfcfld))

	DO J=1,NY
	 DO I=1,NX
	  DO K=1,NSFCFLD
	   SLPOUT(I,J,K)=SLP(I,J,K)
	  ENDDO

	  DO K=1,NZ
	   HTOUT(I,J,K)=HT(I,J,K)
	   UWOUT(I,J,K)=UW(I,J,K)
	   VWOUT(I,J,K)=VW(I,J,K)
	   RHOUT(I,J,K)=RH(I,J,K)
          ENDDO
         ENDDO
	
  	  DO K=1,NSFCFLD
           SLPOUT(145,J,K)=SLP(1,J,K)
          ENDDO

	  DO K=1,NZ
	   HTOUT(145,J,K)=HT(1,J,K)
	   UWOUT(145,J,K)=UW(1,J,K)
	   VWOUT(145,J,K)=VW(1,J,K)
	   RHOUT(145,J,K)=RH(1,J,K)
          ENDDO
         ENDDO

	endif
	
		

	
c
c
c *** Create output file name.
c
      n=index(outdir,' ')-1

Cmp	make ETA_avn to differentiate from files based on eta grib input
	if (nz .eq. 26) then
      model='ETA_avn'
	elseif (nz .eq. 17) then
      model='ETA_rnl'
      nx=145
	endif
      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
      gproj='LL'
	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'

	if (nx .eq. 145) then
	KGDS(8)=0
	KGDS(2)=145	
	endif


        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')

C
C	MODIFY GDS TO REFLECT FLIPPING OF DATA IN N-S SENSE
C

	if (nz .eq. 26) then
	kgds(4)= -kgds(4)
	kgds(7)= -kgds(7)
	endif
        write(6,*) 'writing kgds ', (kgds(I),i=1,14)
        write(14) kgds

        close(14)
	write(6,*) 'past 14 close'
C
C end new stuff

c
c *** Write isobaric upper air data.
c
	write(6,*) 'starting writes to unit1'

	if (nx .ne. 145) then
      write(1) uw
      write(1) vw
      write(1) ht
      write(1) rh
      write(1) pr
      write(1) slp
	else
	WRITE(1) UWOUT
	WRITE(1) VWOUT
	WRITE(1) HTOUT
	WRITE(1) RHOUT
	WRITE(1) PR
	WRITE(1) SLPOUT
	endif

	if (allocated(uwout)) deallocate(uwout)
	if (allocated(vwout)) deallocate(vwout)
	if (allocated(htout)) deallocate(htout)
	if (allocated(rhout)) deallocate(rhout)
	if (allocated(slpout)) deallocate(slpout)
Cmp	
	write(6,*) 'writing these sfc values to the ETA_avn file...'
	do k=1,nsfcfld
	write(6,*) 'field, value ', k, slp(nx/2,ny/2,k)
	enddo
	do k=1,nsfcfld
	write(6,*) 'field, value ', k, slp(3*nx/4,3*ny/4,k)
	enddo

c
      close(1)
c
      return
      end
c
c===============================================================================
c
      subroutine read_degrib(etafile,nxny,nz,pr
     .                      ,ht,tp,rh,uw,vw,slp,atime)
c
      implicit none
c
      integer nxny,nz,i
c
      real pr(nz),ht(nxny,nz),tp(nxny,nz)
     .      ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz)
     .      ,dummy,slp(nxny,12)
c
      integer kgds(200),kpds(50),len,kerr
     .         ,lenpds,lenkgds,nwords,kpdsl
     .         ,j,k,datsav(9)
     .	       ,IRETO,IRET1,JPDS(200),JGDS(200),KNUM

	logical BITMAP(nxny)
c
      character*(*) etafile
      character(LEN=1)::   pds(50)
      character(LEN=12)::   atime
c_______________________________________________________________________________
c
      len=index(etafile//' ',' ')-1
c
c
c *** Fill time stamp (assume all records are for same time).
c

	write(6,*) 'calling get_gds in read_degrib ',etafile
	call get_gds(etafile(1:len),datsav,kpds)
        write(6,*) 'return calling get_gds in read_degrib'

	if (kpds(8) .eq. 100) kpds(8)=0
      write(atime,'(i2.2,i2.2,i2.2,i2.2,i4.4)') 
     .   kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)
c
c *** Fill a missing value flag into first space of each variable.
c
      do k=1,nz
	 ht(1,k)=-99999.
 	 tp(1,k)=-99999.
	 rh(1,k)=-99999.
	 uw(1,k)=-99999.
	 vw(1,k)=-99999.
      enddo

Cmp initialize surface fields to -99999.  so the interp code can handle
Cmp     appropriately

	write(6,*) 'nxny,datsav(1)*datsav(2): ', nxny,datsav(1)*datsav(2)
	do k=1,12
	do j=1,datsav(1)*datsav(2)
	slp(j,k)=-99999.
	enddo
	enddo

Cmp

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

C	write(6,*) 'here2 : calling baopen ', etafile
C	call baopen(11,etafile,IRETO)
	write(6,*) 'here3'
	if (IRETO .ne. 0) write(6,*) 'BAOPEN TROUBLE!!!! ', IRETO

	jpds=-1
	jgds=-1

	jpds(5)=81
	jpds(6)=1
	jpds(7)=0

	call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,1),IRET1)

	write(6,*) 'first GETGB, IRET1= ', IRET1
	if (IRET1 .eq. 0) then
	write(6,*) 'LAND/SEA READ!!!!! '
	write(6,*) (slp(j,1),j=nwords/2,nwords/2+5)
	endif

	jpds(5)=2
	jpds(6)=102
	jpds(7)=0

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,2),IRET1)

	write(6,*) 'PMSL READ!!!!! '
	write(6,*) (slp(j,2),j=nwords/2,nwords/2+8)

	jpds(5)=1
	jpds(6)=1
	jpds(7)=0

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,3),IRET1)

	write(6,*) 'PSFC READ!!!!! '
	write(6,*) (slp(j,3),j=nwords/2,nwords/2+8)

	jpds(5)=7
	jpds(6)=1
	jpds(7)=0

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,4),IRET1)

	write(6,*) 'ZSFC READ!!!!! '
	write(6,*) (slp(j,4),j=nwords/2,nwords/2+8)

	
Cmp     SOIL FIELDS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C

C	AVN stores soil T in 11, ETA in 85
	jpds(5)=11
	jpds(6)=112
	jpds(7)=10


       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,5),IRET1)

	if (IRET1.eq.0) then
	 write(6,*) 'found soil temp at ',jpds(7)
	write(6,*) (slp(j,5),j=nwords/4,nwords/4+8)
	write(6,*) (slp(j,5),j=nwords/2,nwords/2+8)
	endif

	jpds(7)=2600
       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,7),IRET1)
	if (IRET1.eq.0)  then
	write(6,*) 'found soil temp at ',jpds(7)
	write(6,*) (slp(j,7),j=nwords/2,nwords/2+8)
	write(6,*) (slp(j,7),j=nwords/2,nwords/2+8)
	endif

	jpds(7)=10340
	       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,9),IRET1)
	if (IRET1.eq.0)  then
	write(6,*) 'found soil temp at ',jpds(7)
	write(6,*) (slp(j,9),j=nwords/2,nwords/2+5)
	endif

	jpds(7)=25800
       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,11),IRET1)
	if (IRET1.eq.0) then
	write(6,*) 'found soil temp at ',jpds(7)
	write(6,*) (slp(j,11),j=nwords/4,nwords/4+5)
	endif


C	10-200 layer

	jpds(7)=2760
	call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,7),IRET1)
	if (IRET1.eq.0) then
	write(6,*) 'found soil temp at ',jpds(7)
	write(6,*) (slp(j,7),j=nwords/4,nwords/4+8)
	endif

		
		
C       MOISTURE

	jpds(5)=144
	jpds(7)=10

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,6),IRET1)
	if (IRET1.eq.0) then
	 write(6,*) 'found soil wet at ',jpds(7)
        write(6,*) (slp(j,6),j=nwords/4,nwords/4+8)
	endif

        jpds(7)=2600
       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,8),IRET1)
        if (IRET1.eq.0) then
	 write(6,*) 'found soil wet at ',jpds(7)
        write(6,*) (slp(j,8),j=nwords/2,nwords/2+8)
	endif

        jpds(7)=10340
       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,10),IRET1)
        if (IRET1.eq.0) then
	 write(6,*) 'found soil wet at ',jpds(7)
        write(6,*) (slp(j,10),j=nwords/2,nwords/2+5)
	endif

        jpds(7)=25800
       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,12),IRET1)
        if (IRET1.eq.0) then
	 write(6,*) 'found soil wet at ',jpds(7)
        write(6,*) (slp(j,12),j=nwords/2,nwords/2+5)
	endif

C       10-200 layer

        jpds(7)=2760
        call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,8),IRET1)
        if (IRET1.eq.0) then
	 write(6,*) 'found soil wet at ',jpds(7)
        write(6,*) (slp(j,8),j=nwords/4,nwords/4+5)
	endif


	
C 	ISOBARIC DATA

        jpds(6)=100

        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,ht(1,k),IRET1)
        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,tp(1,k),IRET1)
        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,rh(1,k),IRET1)
        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,uw(1,k),IRET1)
        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,vw(1,k),IRET1)
        if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5)

        write(6,*) 'Z,T,Q,U,V ', ht(nxny/2,k),tp(nxny/2,k),
     +  rh(nxny/2,k),uw(nxny/2,k),vw(nxny/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*4 nx,ny,nz,i,j,k,kk
c
      real*4 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===============================================================================

       subroutine get_gds(etafile,gdsinfo,kpds)

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

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

	LOGICAL*1 BITMAP(70000)

        nxny=360*181

        JPDS=-1
        JGDS=-1

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

	write(6,*) 'want to open unit 11 for ', etafile(1:len)
        call baopen(11,etafile(1:len),IRETO)
        write(6,*)'AQUI!!'
        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
        call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,tmp,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)
       write(6,*) gdsinfo

        return
        print *,'GETGDS PROBLEM'
        stop

        end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

                subroutine get_fullgds(etafile,nx,ny,kgds)

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

        integer*4 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

C        write(6,*) 'calling getgb '
C        write(6,*) 'jpds =  ', jpds
C        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

