      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_______________________________________________________________________________
      print*,'========================================================='
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)

       print*,'qqqqqqqqqqqqq',nx,ny,nz
        if (nx .eq.  360 .and. ny .eq. 181) then
        nz=26
        elseif (nx .eq.  720 .and. ny .eq. 361) then !ioannis
        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,16)        !16 was 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*11   atime
      character*7   model

      logical AVN_SOIL2,AVN_SOIL4  ! ioannis
      NAMELIST/FLAGS/AVN_SOIL2,AVN_SOIL4 ! ioannis

       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=16
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
      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

       OPEN (UNIT=47,FILE='../../eta/src/include/soildata.inp'
     $,FORM='FORMATTED')
       READ (47,FLAGS)
       CLOSE (47)
       write(6,*) 'AVN_SOIL2=',AVN_SOIL2
       write(6,*) 'AVN_SOIL4=',AVN_SOIL4

       IF (AVN_SOIL2.and.AVN_SOIL4) THEN
       write(6,*) '!!!!! ERROR in SOIL data'
       write(6,*) '!!!!! AVN_SOIL2 and AVN_SOIL4 are both TRUE'
       write(6,*) '!!!!! ABORT DECODING'
       STOP 777
       ENDIF

       IF ((.NOT.AVN_SOIL2).and.(.NOT.AVN_SOIL4)) THEN
       write(6,*) '!!!!! ERROR in SOIL data'
       write(6,*) '!!!!! AVN_SOIL2 and AVN_SOIL4 are both FALSE'
       write(6,*) '!!!!! ABORT DECODING'
       STOP 778
       ENDIF

       do k=5,nsfcfld

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

cioannis      elseif(k.eq.10.or.k.eq.12.) then
      elseif(k.eq.10.or.k.eq.12.or.k.eq.14.or.k.eq.16) then

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

      do j=1,ny
      do i=1,nx
       IF (AVN_SOIL2) slp(i,j,k)=slp(i,j,8)  ! ioannis
       IF (AVN_SOIL4) slp(i,j,k)=slp(i,j,12) ! ioannis
      enddo
      enddo
      endif
      endif
Cmp
      if (mod(k,2) .eq. 0) then
      do j=1,ny
      do i=1,nx
cioannis      if (slp(i,j,k) .eq. 0) then
      if ((slp(i,j,k).eq.0.0).or.(slp(i,j,k).eq.-99999.)) then ! ioannis
       slp(i,j,k)=0.14
      endif
      enddo
      enddo

      else

      do j=1,ny
        do i=1,nx
cioannis        if (slp(i,j,k) .eq. 0) then
        if ((slp(i,j,k).eq.0.0).or.(slp(i,j,k).eq.-99999.)) then ! ioannis
        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...
cioannis      nsfcfld=12
      nsfcfld=16
      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
      integer nsfcfld ! ioannis
c
      real pr(nz),ht(nxny,nz),tp(nxny,nz)
     .      ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz)
     .      ,dummy,slp(nxny,16)
cioannis     .      ,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=11)::   atime
      integer month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/ !cspir time
c___________________________________________________________________________
c
      len=index(etafile//' ',' ')-1
c
      nsfcfld=16 ! ioannis
c
c *** Fill time stamp (assume all records are for same time).
c

      write(6,*) 'calling get_gds in read_degrib'
      call get_gds(etafile(1:len),datsav,kpds)
      if (kpds(8) .eq. 100) kpds(8)=0
C-----------cspir "all is right with the time" ----------------
C----------- Fix the time according the 18UTC LAPS analysis----
cpk        if ((kpds(14).ne.0) .or. (kpds(14).ne.3)) then
        if (kpds(14).gt.3) then
           print*,'kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)',
     &             kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)
           kpds(11)=kpds(11)+6  !fix the analysis time
cpk           if(kpds(14).gt.03) then !fix the fcst time
            kpds(14)=kpds(14)-6 
cpk           endif
        endif
cpk         if(kpds(10).gt.month(kpds(9))) then
cpk          kpds(10)=kpds(10)-month(kpds(9))
cpk          kpds(9)=kpds(9)+1
cpk         endif
cpk         if(kpds(9).gt.12) then
cpk            kpds(9)=kpds(9)-12
cpk            kpds(8)=kpds(8)+1
cpk         endif
      write(atime,'(i2.2,i2.2,i2.2,i2.2,i3.3)') 
     .   kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)
      print*,'atime: ',atime
C-----------cspir "all is right with the time" ----------------
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)
cioannis      do k=1,12
      do k=1,nsfcfld
      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)
cioannis       real tmp(70000)
       real tmp(261000)

cioannis       LOGICAL*1 BITMAP(70000)
       LOGICAL*1 BITMAP(261000)

cioannis        nxny=360*181
        nxny=720*361

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

