      program degrib2model_driver
c
      implicit none

	integer filemax
	parameter (filemax=30)
c
      integer nx,ny,nz,ntile,l,I,tsave(filemax),  !Degrib grid dimensions
     +	        ICOUNT,JCOUNT,ISAVE,datsav(9),n,tileone(9),iun
     +,ITOT,JTOT,istart,jstart,II,tmin,tmax,J,thold,kpds(200)
c
      character(LEN=255):: degrib_dir,grib_file(filemax),outdir,tname
      character(LEN=255):: ahold
	
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
	do I=1,filemax
        read(5,'(a)') grib_file(I)
	if (grib_file(I) .eq. '9999') then
	 ntile=I-1
	 goto 1103
	endif
	enddo
 1103	continue
      read(5,'(a)') outdir

	write(6,*) 'will process ', ntile, ' tiles'
      l=index(grib_file(1),' ')-1
	write(6,*) 'these are the tile numbers: '

	tmin=99
	tmax=-99
	DO I=1,NTILE
	tname=grib_file(I)
	read (tname(l-1:l),79) tsave(I)
   79	format(I2)
	if (tsave(I).gt.tmax) tmax=tsave(I)
	if (tsave(I).lt.tmin) tmin=tsave(I)
	write(6,*) tsave(I)
	ENDDO

Cmp	sort the grid numbers from smallest to largest (at this
Cmp	point it is important that the SW corner be first in list
Cmp	and NE corner be last in list)

C	write(6,*) 'min value, current first grid ', tmin, tsave(1)
C	write(6,*) 'max value, current last grid ', tmax, tsave(NTILE)

  
C	write(6,*) 'original order ', (tsave(I),I=1,ntile)
	do while(tmin.ne.tsave(1).or.tmax.ne.tsave(NTILE))
	write(6,*) 'sorting the tile numbers'

	do J=1,ntile-1
	if (tsave(J+1).lt.tsave(J)) then
	thold=tsave(J)
	ahold=grib_file(J)
	tsave(J)=tsave(J+1)
	grib_file(J)=grib_file(J+1)
	tsave(J+1)=thold
	grib_file(J+1)=ahold
	endif
	enddo
C	write(6,*) 'first and last grids now ', tsave(1),tsave(ntile)
	write(6,*) 'new order ', (tsave(I),I=1,ntile)

	enddo

c
CCCCCCCCCCCCCCCCCCCCC
       nz=39
!      nz=38
Cmp
CCCCCCCCCCCCCCC
c

	ITOT=0.
	JTOT=0.

	DO I=1,NTILE

	tname=grib_file(I)
	n=index(tname,' ')-1

	iun=11+I
        call get_gds(tname(1:n),datsav,kpds,iun)
C	write(6,*) 'datsav on tile: ', I, ' is ',(datsav(II),II=1,5)

	if (I .eq. 1) then
	  do II=1,9
	    tileone(II)=datsav(II)
	  enddo
	endif

	call get_index(tileone,istart,jstart,datsav)

	ITOT=max(istart-1,ITOT)
	JTOT=max(jstart-1,JTOT)

	if (I .eq. ntile) write(6,*) 'computed NX,NY= ', ITOT+datsav(1),
     + JTOT+datsav(2)
	nx=ITOT+datsav(1)
	ny=JTOT+datsav(2)

	ENDDO	

       call degrib2model(grib_file,outdir,nx,ny,nz,ntile)

      end
c
c===============================================================================
c
      subroutine degrib2model(file_in,outdir,nx,ny,nz,ntile)
c
      implicit none
c
      real    cp,kappa
      parameter(cp=1004.,kappa=287./1004.)
c
      integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld,ntile
     .     ,iyear,imonth,iday,ifcsthr
     .     ,tile,ival,jval,istart,jstart,icount,jcount
     .	   ,tilenum,tmpval,datsav(9),istep,jstep,tileone(9)
     .     ,imax,jmax,III,IRETO,kpds(200),iun,kgds(200)
c
      real ht(nx,ny,nz),htin(59,47,nz) !Isobaric heights (m)
     .  ,tp(nx,ny,nz),tpin(59,47,nz) !Isobaric temps (K)
     .  ,th(nx,ny,nz),thin(59,47,nz) !Isobaric theta (K)
     .  ,uw(nx,ny,nz),uwin(59,47,nz) !Isobaric u-wind (m/s)
     .  ,vw(nx,ny,nz),vwin(59,47,nz) !Isobaric v-wind (m/s)
     .  ,rh(nx,ny,nz),rhin(59,47,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
     . ,slp(nx,ny,12),slpin(59,47,12),llcorner(2,30)
     . ,tmp(nx*ny,nz)

	real urlat,urlon
	real tbar,old,fact,biassum,z1000t,z975t
c
Cmp      character*255 degrib_dir,degrib_file,outdir,outfile,file_in(ntile)
      character(LEN=255):: degrib_dir,degrib_file,outdir,outfile
      character(LEN=255):: gdsfile,file_in(30)
      character(LEN=2)::   gproj
      character(LEN=11)::   atime
      character(LEN=8)::   model
c
      common /estab/esat(15000:45000),es(15000:45000)

	common /maxval/ imax,jmax

Cmp Loop over whole subroutine

        do III=1,ntile
        degrib_file=file_in(III)
        tilenum=III

c_______________________________________________________________________________
c
c *** Fill pressure levels.
c
Cmp
	nsfcfld=12

Cmp	pressure levels are hardwired here

      do k=1,nz
         pr(k)=1025.-float(k*25)
      enddo
c 
c *** Read in degrib data.
c
       l=index(degrib_file,' ')-1

	iun=21+III
        call get_gds(degrib_file(1:l),datsav,kpds,iun)

        call read_degrib(degrib_file(1:l)
     .  ,datsav(1)*datsav(2),nz,pr,htin,tpin,rhin,uwin,vwin,slpin,atime)


	if (tilenum .eq. 1) then
	do I=1,9
	 tileone(I)=datsav(I)
	enddo
	endif

cmp	Want to do computations relative to the point tileone(3),tileone(4)
	
	llcorner(1,tilenum)=datsav(3)
	if (datsav(4) .ge. 180000) datsav(4)=datsav(4)-360000
	llcorner(2,tilenum)=datsav(4)
	
Cmp	The istep, jstep, and index values should be obtained from
Cmp	the datsav data in comparison with the tile 1 datsav data.

	call get_index(tileone,istart,jstart,datsav)

	write(6,*) 'writing from I,J= ', istart,jstart
C	write(6,*) 'writing out for J= ', jstart,jstart+datsav(2)-1
	do K=1,nz
	do I=istart,istart+datsav(1)-1
	do J=jstart,jstart+datsav(2)-1
	ival=(I-istart)+1
	jval=(J-jstart)+1
	if (ival .gt. 59 .or. ival .lt. 1 .or. jval .gt. 47)
     +    write(6,*) 'ival,jval= ', ival,jval
	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

	do K=1,nsfcfld
	do I=istart,istart+datsav(1)-1
        do J=jstart,jstart+datsav(2)-1
	ival=(I-istart)+1
        jval=(J-jstart)+1
	slp(I,J,K)=slpin(ival,jval,K)
	enddo
	enddo
	enddo

C	write(6,*) 'check some vals along a seam', istart-6, istart+2
	do J=jstart,jstart-8,-1
C	write(6,298) (rh(I,J,13),I=istart-6,istart+2)
	enddo
C	write(6,*) 'soil temps'
	do J=jstart,jstart-8,-1
C	write(6,298) (slp(I,J,5),I=istart-6,istart+2)
	enddo
  298	format(9(f7.2,1x))


C	write(6,*) 'SLP values '
	do K=1,12
C	write(6,*) slp(istart+datsav(1)/2,jstart+datsav(2)/2,K)
	enddo

Cmp
Cmp	Handle missing surface data....
Cmp
	do k=1,12
	if (slp(1,1,k) .eq. -99999.) then

	write(6,*) 'SURFACE DATA MISSING... ', K

      if (k.eq.9.or.k.eq.11.) then
C	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
C	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
	enddo

c *** Convert 3d temp to theta.
c *** Compute Exner function.
c *** Convert 3d rh to mr.
c
Cmp
Cmp	Only do this stuff when the grid is completely filled!
Cmp
	if (tilenum .eq. ntile) then

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

Cnew
C***
C***    Change 1000 hPa heights to match lowest temps
C***
        biassum=0.

        do J=1,NY
        do I=1,NX

        fact=287.*(pr(2)-pr(3))/(9.81*pr(2))
        z975t=(ht(i,j,3)-ht(i,j,2))/fact
        fact=287.*(pr(1)-pr(2))/(9.81*pr(1))
        z1000t=(ht(i,j,2)-ht(i,j,1))/fact

C
Cthis criteria empirically derived.  Was found that where 950-975 thick >
Cabout 7 m greater than 975-1000 thick, lowest-lev temperatures became hosed.
C
        if (z975t-z1000t .gt. 2.5) then
C       tbar=(1./4.)*(tp(i,j,1)+tp(i,j,2)+tp(i,j,3)+tp(i,j,4))
        tbar=amax1(tp(i,j,1),tp(i,j,2),tp(i,j,3),tp(i,j,4))
        ht(i,j,1)=ht(i,j,2)-fact*tbar
        biassum=biassum+1
        endif

        enddo
        enddo

	write(6,*) 'modified Z1000 at ', biassum, ' points'

Cendnew


      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
c
      print *,'ua :',ht(nx/2,ny/2,nz/2),th(nx/2,ny/2,nz/2),
     .	rh(nx/2,ny/2,nz/2)
     .       ,uw(nx/2,ny/2,nz/2),vw(nx/2,ny/2,nz/2)
c
c *** Create output file name.
c
	write(6,*) 'WRITING THE DATA TO FILE'

      n=index(outdir,' ')-1
      model='ETA_tile'
      outfile=outdir(1:n)//'/'//atime//'.'//model
      n=index(outfile,' ')-1
      open(1,file=outfile(1:n),status='unknown',form='unformatted')
c
c *** Write header stuff.
c
      ip=1
      jp=1
      nsfcfld=12
      gproj='LC'
!      write(1) nx,ny,nz,nx,ny,ip,jp,nsfcfld,gproj
	write(1) nz

	lat1=datsav(8)/1000.
	lat2=datsav(9)/1000.
	lon0=datsav(5)/1000.

Cmp
Cmp	The sw(1) and sw(2) values are passed from the degribbing, use the
Cmp	lambert subroutine to calculate the ne values (assumes that the
Cmp	most NE'ly tile was processed last and provides the GDS info in
Cmp	datsav)
Cmp
Cmp	Could this obtaining of NE corner point be made more general so
Cmp	as to prevent needing to process that grid last?

	write(6,*) 'gds based sw values ', llcorner(1,1)/1000.,
     + llcorner(2,1)/1000.
	call lambert (urlat,urlon,datsav)
	write(6,*) 'gds based ne corner ', urlat,urlon

	sw(1)=llcorner(1,1)/1000.
	sw(2)=llcorner(2,1)/1000.
	ne(1)=urlat
	ne(2)=urlon

        n=index(degrib_file,' ')-1

        call get_fullgds(degrib_file(1:n),datsav(1),datsav(2),kgds)

C	overwrite grid dims and sw corner coordinates

	kgds(2)=nx
	kgds(3)=ny
	kgds(4)=sw(1)*1000.
	kgds(5)=sw(2)*1000.

C new stuff for GDS file
C
        n=index(outdir,' ')-1
        gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model

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

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

        close(14)
C
C end new stuff


!     write(1) nx,ny,nz,lat1,lat2,lon0,sw,ne
c
c *** Write isobaric upper air data.
c

Cmp	Write out every 10th grid point in both directions?

Cnot	goto 1076
	write(6,*) 'level 20', nx,ny

	write(6,*) 'mixr * 1000.'
	do J=ny,1,-(ny/15)
	write(6,744) (rh(I,J,20)*1000.,I=1,nx,nx/9)
	enddo

	write(6,*) 'HT'
	do J=ny,1,-(ny/15)
	write(6,744) (ht(I,J,20),I=1,nx,nx/9)
	enddo

	write(6,*) 'UW'
	do J=ny,1,-(ny/15)
        write(6,744) (uw(I,J,20),I=1,nx,nx/9)
        enddo

	write(6,*) 'VW'
	do J=ny,1,-(ny/15)
        write(6,744) (vw(I,J,20),I=1,nx,nx/9)
        enddo

	write(6,*) 'TH'
	do J=ny,1,-(ny/15)
        write(6,744) (th(I,J,20),I=1,nx,nx/9)
        enddo

 1076	continue

  743	format(30(e9.3,1x))
  744   format(30(f7.1,1x))

      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

        do k=1,nsfcfld
        write(6,*) 'field, value ', k, slp(nx/2,ny/2,k)
        enddo
c
      close(1)
	endif
c
        enddo

	write(6,*) 'leaving degrib2model'
       return
       end
c
c===============================================================================
c
C      subroutine read_degrib(etafile,nxny,nz,pr
C     .                      ,ht,tp,rh,uw,vw,slp,atime)
      subroutine read_degrib(etafile,nxny,nz,pr
     .       ,htout,tpout,rhout,uwout,vwout,slpout,atime)
c
      implicit none
c
      integer nxny,nz,i,imax,jmax
c
      real pr(nz),ht(nxny,nz),tp(nxny,nz)
     .      ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz)
     .      ,dummy,slp(nxny,12)

      real htout(imax,jmax,nz),tpout(imax,jmax,nz),rhout(imax,jmax,nz),
     .     uwout(imax,jmax,nz),vwout(imax,jmax,nz),slpout(imax,jmax,12)

	real crot(nxny),srot(nxny),rmagb,rmagaft,ubef,vbef,urlat,
     .      urlon

c
      integer kgds(200),kpds(200),len,kerr
     .         ,lenpds,lenkgds,nwords,kpdsl
     .         ,j,k,datsav(9)
     .	       ,jpds(200),jgds(200),IRETO,KNUM,IRET1,iun
c
      character*(*) etafile
      character(LEN=11)::   atime

	logical BITMAP(nxny)

	common /maxval/imax,jmax
c_______________________________________________________________________________
c
       len=index(etafile//' ',' ')-1
c
C	write(6,*) 'inside read_degrib, dimension with ', nxny,nz

	iun=31
        call get_gds(etafile(1:len),datsav,kpds,iun)

        if (kpds(8) .ge. 100) kpds(8)=kpds(8)-100

      write(atime,'(i2.2,i2.2,i2.2,i2.2,i3.3)')
     .   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 

	do k=1,12
	do j=1,nxny
	slp(j,k)=-99999.
	enddo
	enddo

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

1     continue

C  add something to read in surface fields in here
C
        call baopen(11,etafile,IRETO)
        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)

C        if (IRET1 .eq. 0) write(6,*) 'LAND/SEA READ!!!!! '
C        write(6,*) (slp(j,1),j=nwords/2,nwords/2+3)

        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)

C	if (IRET1 .eq. 0) write(6,*) 'PMSL READ!!!!! '
C        write(6,*) (slp(j,2),j=nwords/2,nwords/2+3)

        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)

C        if (IRET1 .eq. 0) write(6,*) 'PSFC READ!!!!! '
C        write(6,*) (slp(j,3),j=nwords/2,nwords/2+3)

        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)

c       	if (IRET1 .eq. 0)  write(6,*) 'ZSFC READ!!!!! '
C        write(6,*) (slp(j,4),j=nwords/2,nwords/2+3)

C
Cmp     SOIL FIELDS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
	
C       if (mod(kpds(7)-i,256).eq.0) then
Cmp     write(6,*) 'kpds(5)= ', kpds(5)
Cmp     write(6,*) 'lower is ', i
Cmp     write(6,*) 'upper is ', (kpds(7)-i)/256.
C       endif

        jpds(5)=85
        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
c        write(6,*) 'found soil temp over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,5),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil temp over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,7),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil temp over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,9),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil temp over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,11),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil wet over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,6),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil wet over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,8),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil wet over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,10),j=nwords/2,nwords/2+3)
	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
C        write(6,*) 'found soil wet over ',jpds(7), 'layer!!'
C        write(6,*) (slp(j,12),j=nwords/2,nwords/2+3)
	endif

C       jpds(5)=7
        jpds(6)=100

Cmp     doing this in proper order??????
Corig   do k=1,nz
        do k=nz,1,-1
        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)

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

        enddo
c
c *** Normal finish.
c
1000  continue

Cmp     rotate the winds at this point
C	write(6,*) 'rotating a lambert projection'
        call rotate_lcc(kgds,crot,srot,nxny)

	do K=1,nz
        do I=1,nxny
        rmagb=(uw(I,K)**2. + vw(I,K)**2.)**(0.5)
        ubef=uw(I,K)
        vbef=vw(I,K)

        uw(I,K)=crot(I)*ubef+srot(I)*vbef
        vw(I,K)=crot(I)*vbef-srot(I)*ubef
	rmagaft=(uw(I,K)**2. + vw(I,K)**2.)**(0.5)

        if (abs(rmagaft-rmagb).gt.3.) then
        write(6,*) 'MAG, I,K,old,new==> ',I,K,rmagb,rmagaft
        write(6,*) 'original components..', ubef,vbef
        write(6,*) 'new components..', uw(I,k),vw(I,K)
        write(6,*) 'rotation cosines ', crot(I),srot(I)
        write(6,*) '.................................'
        endif
        ENDDO
        ENDDO

Cmp resort the data into the output arrays at this point

	call resort_1d(ht,htout,nz,datsav)
	call resort_1d(uw,uwout,nz,datsav)
	call resort_1d(vw,vwout,nz,datsav)
	call resort_1d(tp,tpout,nz,datsav)
	call resort_1d(rh,rhout,nz,datsav)
	call resort_1d(slp,slpout,12,datsav)

      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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	        subroutine lambert(gdslatur,gdslonur,datsav)

C
C       Subroutine written 9 March 1999 by M. Pyle to support tiled 221 input.
C       Code adapted from GEMPAK routine gblamb.c.  Whole purpose is to get the
C       UR corner lat and lon for use in workstation Eta.

        integer latin1,latin2,nx,ny,la1,lo1,lov
        integer dx,dy,datsav(9)

        real(kind=8) earth_rad, const_cone, xll,yll,xur,yur,lat1,lon1,loncnt,
     +  angle1,angle2,x1,x2,y1,y2,alpha,rtemp
	real gdslatur,gdslonur
     +  ,gdslatll,gdslonll

        parameter(rpi=3.141592654)
        parameter(d2r=rpi/180.)
        parameter(r2d=180./rpi)
        parameter(radius=6370000.)

Ctest        common /gdsinfo/ datsav(9)

        latin1=datsav(8)
        latin2=datsav(9)
        la1=datsav(3)
        lo1=datsav(4)
        dx=datsav(6)
        dy=datsav(7)
        nx=datsav(1)
        ny=datsav(2)
        lov=datsav(5)

C        write(6,*) 'values in lambert '
C        write(6,*) latin1,latin2,la1,lo1,dx,dy,nx,ny,lov

        lat1= (la1/1000.0)*d2r
        if (lo1 .eq.  180000) lo1=lo1-360000
        if (lo1 .lt. -180000) lo1=lo1+360000
        lon1= (lo1/1000.0)*d2r

Cmp     now have LL corner in radians, W is negative

        if (lov .eq.  180000) lov=lov-360000
        if (lov .lt. -180000) lov=lov+360000
        loncnt= (lov/1000.0)*d2r

        angle1= (rpi/2.) - ( abs(latin1/1000.0) * d2r )
        angle2= (rpi/2.) - ( abs(latin2/1000.0) * d2r )

        if (latin1 .eq. latin2) then
        const_cone=cos(angle1)
        else
        const_cone= ( log ( sin(angle2) ) - log ( sin ( angle1 ) ) )/
     +  ( log ( tan ( angle2/2 ) ) - log ( tan ( angle1/2 ) ) )
        endif

C        write(6,*) 'const_cone= ', const_cone

        earth_rad=radius/const_cone

cmp     assuming NH

        x1 = earth_rad * tan( (rpi/2.-lat1) / 2 )**(const_cone)*
     +  sin (const_cone * ( lon1 - loncnt ) )
        y1 = -earth_rad * tan( (rpi/2.-lat1) / 2 )**(const_cone)*
     +  cos (const_cone * ( lon1 - loncnt ) )

	 alpha= (tan(angle1 / 2 )**const_cone)/sin (angle1)

        x2=x1 + ( nx - 1 ) * alpha * dx
        y2=y1 + ( ny - 1 ) * alpha * dy

        xll=min(x1,x2)
        xur=max(x1,x2)
        yll=min(y1,y2)
        yur=max(y1,y2)

	xlltmp=abs(xll)
	ylltmp=abs(yll)
        gdslatll= ( rpi/2. - 2 *
     + atan ( ( (xlltmp**2.+ylltmp**2.)**(0.5)/earth_rad)**
     + (1/const_cone) ) ) * r2d
	

        rtemp= atan2 ( xll, -yll ) * ( 1 / const_cone ) + loncnt

        if ( rtemp .gt. rpi ) then
        gdslonll = ( rtemp - 2.*rpi ) * r2d
        else if ( rtemp .lt. -rpi ) then
        gdslonll = ( rtemp + 2.*rpi ) * r2d
        else
        gdslonll = rtemp * r2d
        endif

	xurtmp=abs(xur)
	yurtmp=abs(yur)

        gdslatur= ( rpi/2. - 2 *
     + atan ( ( (xurtmp**2.+yurtmp**2.)**(0.5)/earth_rad)**
     + (1/const_cone) ) ) * r2d

        rtemp= atan2 ( xur, -yur ) * ( 1 / const_cone ) + loncnt
        if ( rtemp .gt. rpi ) then
        gdslonur = ( rtemp - 2.*rpi ) * r2d
        else if ( rtemp .lt. -rpi ) then
        gdslonur = ( rtemp + 2.*rpi ) * r2d
        else
        gdslonur = rtemp * r2d
        endif

C        write(6,*) 'output==> '
C        write(6,*) 'LL points ', gdslatll,gdslonll
C        write(6,*) 'UR points ', gdslatur,gdslonur

	return

        END

C***************************************************************
C
	subroutine resort_1d(arrayin,arrayout,nz,datsav)

	integer imax,jmax,nz,datsav(9)

	common /maxval/imax,jmax

	real arrayout(imax,jmax,nz)
	real arrayin(datsav(1)*datsav(2),nz)

	arrayout=-9999.

	do 90 K=1,nz

	do J=1,datsav(2)
	do I=1,datsav(1)
	indval=(J-1)*datsav(1)+I
	arrayout(I,J,K)=arrayin(indval,K)
	enddo
	enddo

   90	continue

	return
	end
C*************************************************************

        subroutine rotate_lcc(kgds,crot,srot,npts)

Cmp     stolen/adapted from iplib code gdswiz03
C
C SUBPROGRAM:  GDSWIZ03   GDS WIZARD FOR LAMBERT CONFORMAL CONICAL
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION
C           (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63)
C           AND RETURNS ONE OF THE FOLLOWING:
C             (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
C             (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
C           FOR LAMBERT CONFORMAL CONICAL PROJECTIONS.
C           IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
C           BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
C           OUTPUT ELEMENTS ARE SET TO FILL VALUES.
C           THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.

C       LAMBERT CONFORMAL GRIDS
C          (2)   - NX NR POINTS ALONG X-AXIS
C          (3)   - NY NR POINTS ALONG Y-AXIS
C          (4)   - LA1 LAT OF ORIGIN (LOWER LEFT)
C          (5)   - LO1 LON OF ORIGIN (LOWER LEFT)
C          (6)   - RESOLUTION (RIGHT ADJ COPY OF OCTET 17)
C          (7)   - LOV - ORIENTATION OF GRID
C          (8)   - DX - X-DIR INCREMENT
C          (9)   - DY - Y-DIR INCREMENT
C          (10)  - PROJECTION CENTER FLAG
C          (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
C          (12)  - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER
C          (13)  - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER


C        parameter (NPTS=datsav(1)*datsav(2))
      INTEGER KGDS(200)
      REAL RLON(NPTS),RLAT(NPTS)
      REAL CROT(NPTS),SROT(NPTS)
        real DLON,AN
        real  DE,DR
        REAL PI,DPR
      PARAMETER(RERTH=6.3712E6)
      PARAMETER(PI=3.14159265,DPR=180./PI)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

	FILL=-999.
        LROT=1
        IROT=1
        IM=KGDS(2)
        JM=KGDS(3)
        RLAT1=KGDS(4)*1.E-3
        RLON1=KGDS(5)*1.E-3
        IROT=MOD(KGDS(6)/8,2)
        ORIENT=KGDS(7)*1.E-3
        DX=KGDS(8)
        DY=KGDS(9)
        IPROJ=MOD(KGDS(10)/128,2)
        ISCAN=MOD(KGDS(11)/128,2)
        JSCAN=MOD(KGDS(11)/64,2)
        NSCAN=MOD(KGDS(11)/32,2)
        RLATI1=KGDS(12)*1.E-3
        RLATI2=KGDS(13)*1.E-3
        H=(-1.)**IPROJ
        HI=(-1.)**ISCAN
        HJ=(-1.)**(1-JSCAN)
        DXS=DX*HI
        DYS=DY*HJ

        IF(RLATI1.EQ.RLATI2) THEN
          AN=SIN(H*RLATI1/DPR)
        ELSE
          AN=LOG(COS(RLATI1/DPR)/COS(RLATI2/DPR))/
     &       LOG(TAN((H*RLATI1+90)/2/DPR)/TAN((H*RLATI2+90)/2/DPR))
        ENDIF        
	DE=RERTH*COS(RLATI1/DPR)*TAN((H*RLATI1+90)/2/DPR)**AN/AN
        IF(H*RLAT1.EQ.90) THEN
          XP=1
          YP=1
        ELSE
          DR=DE/TAN((H*RLAT1+90)/2/DPR)**AN
          DLON1=MOD(RLON1-ORIENT+180+3600,360.)-180
C       atmp=RLON1-ORIENT+180.+3600.
C       DLON1= atmp - INT (atmp/360.)*360. - 180.
          XP=1-H*SIN(AN*DLON1/DPR)*DR/DXS
          YP=1+COS(AN*DLON1/DPR)*DR/DYS
        ENDIF
        ANTR=1/(2*AN)
        DE2=DE**2
        XMIN=0
        XMAX=IM+1
        YMIN=0
        YMAX=JM+1
        NRET=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
C       XP=1
C       YP=1
        DO N=1,IM*JM
        J=INT((N-1)/IM)+1
        I=N-(J-1)*IM
            IF(I.GE.XMIN.AND.I.LE.XMAX.AND.
     &         J.GE.YMIN.AND.J.LE.YMAX) THEN
              DI=(I-XP)*DXS
              DJ=(J-YP)*DYS
              DR2=DI**2+DJ**2
              IF(DR2.LT.DE2*1.E-6) THEN
                RLON(N)=0.
                RLAT(N)=H*90.
              ELSE
                RLON(N)=MOD(ORIENT+H/AN*DPR*ATAN2(DI,-DJ)+3600,360.)
C       atmp=ORIENT+H/AN*DPR*ATAN2(DI,-DJ)+3600
C       RLON(N)= atmp - INT(atmp/360.) * 360.
                RLAT(N)=H*(2*DPR*ATAN((DE2/DR2)**ANTR)-90)
              ENDIF
              NRET=NRET+1
              IF(LROT.EQ.1) THEN
                IF(IROT.EQ.1) THEN
C       atmp=RLON(N)-ORIENT+180.+3600.
C       DLON=atmp - INT(atmp/360.)*360.-180.
                  DLON=MOD(RLON(N)-ORIENT+180+3600,360.)-180
                  CROT(N)=H*COS(AN*DLON/DPR)
                  SROT(N)=SIN(AN*DLON/DPR)
                ELSE
                  CROT(N)=1
                 SROT(N)=0
                ENDIF
              ENDIF
            ELSE
              RLON(N)=FILL
              RLAT(N)=FILL
            ENDIF
          ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        RETURN
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	BLOCK DATA DIMS
	common /maxval/ imax,jmax
	data imax/59/
	data jmax/47/
	END BLOCK DATA DIMS

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c===============================================================================

       subroutine get_gds(etafile,gdsinfo,kpds,iun)

        character(LEN=255):: etafile

        integer kgds(200),kpds(200),len
     .         ,nwords
     .         ,j,k,gdsinfo(9),iun,II
     .         ,IRETO,JGDS(200),JPDS(200)
        real tmp(5000)
        logical bitmap(5000)

        nxny=5000
        JPDS=-1
        JGDS=-1
        len=index(etafile//' ',' ')-1

C        write(6,*) 'calling baopen for iun,file: ', iun,etafile(1:len)
        call baopen(iun,etafile(1:len),IRETO)

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

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

        call getgb(iun,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,tmp,IRET1)

        if (IRET1 .ne. 0) then
        write(6,*) 'back from getgb in get_gds: ', IRET1
        endif

       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,*) 'gds in get_gds: ', (gdsinfo(II),II=1,4)

        call baclose(iun,IRET)

        if (IRET .ne. 0) write(6,*) 'bad baclose in get_gds'

        return
        end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


        subroutine get_index(tileone,outi,outj,datsav)

C
C       Subroutine written 11 March 1999 by M. Pyle to support tiled 221 input.
C       Code adapted from GEMPAK routine gblamb.c.
C       Purpose is to do a K. Brill worthy piece of tiling code that will
C       figure out where to put the data in the output data set (get the
C       I and J indices)

        integer latin1,latin2,nx,ny,la1,lo1,lov,tonelat,tonelon
        integer dx,dy,datsav(9),outi,outj,tileone(9)

        real(kind=8):: earth_rad, const_cone, xll,yll,xur,yur,lat1,
     +  angle1,angle2,x1,x2,y1,y2,alpha,rtemp,latfst,lonfst,lon1,loncnt
        real gdslatur,gdslonur
     +  ,gdslatll,gdslonll

        parameter(rpi=3.141592654)
        parameter(d2r=rpi/180.)
        parameter(r2d=180./rpi)
        parameter(radius=6370000.)


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       latin values are the true "cutting" latitudes of the projection
C
C       la1, lo1 are the lat/lon of the LL point of the GRIB data
C
C       dx,dy is a measure of grid spacing
C
C       nx,ny are the I and J dimensions of the GRIB data
C
C       lov is the center longitude of the projection
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc

        latin1=datsav(8)
        latin2=datsav(9)
        la1=datsav(3)
        lo1=datsav(4)
        tonelat=tileone(3)
        tonelon=tileone(4)
        dx=datsav(6)
        dy=datsav(7)
        nx=datsav(1)
        ny=datsav(2)
        lov=datsav(5)

        latfst=(tonelat/1000.0)*d2r
        lat1= (la1/1000.0)*d2r
        if (lo1 .ge.  180000) lo1=lo1-360000
        if (lo1 .lt. -180000) lo1=lo1+360000
        if (tonelon .ge.  180000) tonelon=tonelon-360000
        if (tonelon .lt. -180000) tonelon=tonelon+360000
        lon1= (lo1/1000.0)*d2r
        lonfst= (tonelon/1000.0)*d2r

C       write(6,277) lat1*r2d,lon1*r2d
C       write(6,278) latfst*r2d,lonfst*r2d
  277   format( 'comparing the point ',2(f8.3,1x))
  278   format( 'with tile one   ..> ',2(f8.3,1x))

Cmp     now have LL corner in radians, W is negative

        if (lov .ge.  180000) lov=lov-360000
        if (lov .lt. -180000) lov=lov+360000
        loncnt= (lov/1000.0)*d2r

        angle1= (rpi/2.) - ( abs(latin1/1000.0) * d2r )
        angle2= (rpi/2.) - ( abs(latin2/1000.0) * d2r )

        if (latin1 .eq. latin2) then
        const_cone=cos(angle1)
        else
        const_cone= ( log ( sin(angle2) ) - log ( sin ( angle1 ) ) )/
     +  ( log ( tan ( angle2/2 ) ) - log ( tan ( angle1/2 ) ) )
        endif

C        write(6,*) 'const_cone= ', const_cone

        earth_rad=radius/const_cone

cmp     assuming NH

        x1 = earth_rad * tan( (rpi/2.-latfst) / 2 )**(const_cone)*
     +  sin(const_cone * ( lonfst - loncnt ) )

        x2 = earth_rad * tan( (rpi/2. - lat1) / 2 )**(const_cone)*
     +  sin(const_cone * ( lon1  -  loncnt ) )

        y1=- earth_rad * tan( (rpi/2.-latfst) / 2 )**(const_cone)*
     +  cos(const_cone * ( lonfst - loncnt ) )

        y2 =-earth_rad * tan( (rpi/2. - lat1) / 2 )**(const_cone)*
     +  cos(const_cone * ( lon1  -  loncnt ) )

         alpha= (tan(angle1 / 2 )**const_cone)/sin (angle1)

C       write(6,275) x1,x2
C       write(6,276) y1,y2
  275   format('x vals ', 2(f12.2,1x))
  276   format('y vals ', 2(f12.2,1x))

        xll=min(x1,x2)
        xur=max(x1,x2)
        yll=min(y1,y2)
        yur=max(y1,y2)

        if ( (xur - xll) .lt. 100.) then
Cmp     write(6,*) 'same column as first tile....set outi=1'
        outi=1
        else
        outi= (xur-xll)/(alpha*dx) + 2
        endif

        if ( (yur - yll ) .lt. 100.) then
Cmp     write(6,*) 'same row as first tile....set outj=1'
        outj=1
        else
        outj= (yur-yll)/(alpha*dy) + 2
        endif

C       write(6,*) 'out I and out J= ', outi,outj

        return

        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                subroutine get_fullgds(etafile,nx,ny,kgds)

        character*(*) etafile

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

        nxny=nx*ny

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

        jpds=-1
        jgds=-1

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

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

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

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