      program degrib2model_driver
c
      implicit none
c
      integer nx,ny,nz          !Degrib grid dimensions
c
      character(LEN=255):: degrib_dir,degrib_file,outdir,gribfile
c
      real esat,es
      integer datsav,n,l,kpds(200)
c
      common /estab/esat(15000:45000),es(15000:45000)
      common /gdsinfo/datsav(11)
c_______________________________________________________________________________
c
c *** Initialize tables.
c
      call es_ini
c
c *** Read file names.
c
        read(5,'(a)') gribfile
        read(5,'(a)') outdir

      n=index(gribfile,' ')-1
      write(6,*) 'GRIBFILE is ', gribfile(1:n)
      CALL GET_GDS(gribfile(1:n),datsav,kpds)

c
C     X and Y dimensions now obtained from the degribbed GRIB file
c
       nx=datsav(1)
       ny=datsav(2)
       nz=19
c
       write(6,*) 'calling degrib2model ', nx,ny,nz
       call degrib2model(gribfile,outdir,nx,ny,nz)
c
      end
c
c===============================================================================
c
      subroutine degrib2model(gribfile,outdir,nx,ny,nz)
c
      implicit none
c
!Niv
      logical contAss   !Se for falso, nao ha assimilacao...
!Niv
      real    cp,kappa
      parameter(cp=1004.,kappa=287./1004.)
c
      integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld
     .         ,iyear,imonth,iday,ifcsthr
     .	       ,datsav,ii,jj
      integer kgds(200)

	real crot(nx*ny),srot(nx*ny),tbar,old,fact,biassum
	real z1000t,z975t
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)
     .      ,etalevs(19)
     .      ,xe,mrsat,esat,es
     .	    ,slp(nx,ny,12) 	  !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL
				  !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC
c
!Niv Declaracao das variaveis utilizadas na assimilacao
!
       integer dia,hora,minuto
       integer istatAss1,istatAss2,istatAss3,istatAss4
       integer kmObs,nObs
       parameter(kmObs=100)
       character horAss*2   !identifica arquivo que sera modificado
       character dadosAssim*53,nomeAer(kmObs)*4
       character               aero*4
       integer ind(nx,ny)
       real latGr(ny),lonGr(nx)
       real deltay   ,deltax
       real altiObs(kmObs),latiObs(kmObs),longObs(kmObs)
       real alti,lati,longi
       real pnmmObs(kmObs),tempObs(kmObs)
       real uObs(kmObs),vObs(kmObs)
       real dir,dirRad,NPI,veloc
       real dif
!
!Niv
      character(LEN=255):: degrib_dir,degrib_file
      character(LEN=255):: outdir,outfile,gribfile
      character(LEN=255):: gdsfile
      character(LEN=2)  :: gproj
      character(LEN=12) ::   atime
      character(LEN=6)  ::  model

	real urlat,urlon
c
       data ETALEVS/1000., 925., 900., 850., 800.,
     +               750., 700., 650., 600., 550.,
     +               500., 450., 400., 350., 300,
     +               250., 200., 150., 100./

      common /estab/esat(15000:45000),es(15000:45000)
      common /gdsinfo/datsav(11)
c_______________________________________________________________________________
c
c *** Fill pressure levels.
c
Cmp
	nsfcfld=12

Cmp     pressure levels are hardwired here

      do k=1,nz
         pr(k)=ETALEVS(k)
	 write(6,*) 'PRESS ',pr(k)
      enddo

c
c *** Read in degrib data.
c
      n=index(gribfile,' ')-1
      horAss(1:1)=gribfile(n-5:n-5)
      horAss(2:2)=gribfile(n-4:n-4)
	write(6,*) 'calling read_degrib'


      call read_degrib(gribfile(1:n)
     .                ,nx*ny,nz,pr,ht,tp,rh,uw,vw,slp,atime)

	write(6,*) 'BACK FROM read_degrib'
      write(*,*) 'horAss ',horAss(1:2)

      contAss=.true.
!      contAss=.false.
      if(horAss(1:2).eq.'00'.and.(contAss)) then

      write(99,*) 'gribfile ',gribfile(1:n)
      write(99,*)'horAss ',horAss(1:2)
!Niv
!     21 out 2004
!     Esquema de Assimilacao - Vou inserir a partir deste trecho.
!     Assimilando: Metar
!                  kmObs  - quantidade de observacoes
!     localidades: Os aeroportos observados estao relacionados
!                  na pasta: /home2/nivaldo/worketa_all/data/metar/saida
!                  arquivo: locais.txt
!                  nomAer(kmObs)
!                  longObs(kmObs)
!                  latiObs(kmObs)
!                  altiObs(kmObs)
!     Os dados estao na pasta: /home2/nivaldo/worketa_all/data/metar/2004/mm/dd/hh.txt
!                  mm - mes
!                  dd - dia
!                  hh.txt - hora da observacao
!                  00.txt - corresponde a zero hora
!     Variaveis:
!                  PNMM   - pnmmObs(kmObs)
!                  TEMP   - tempObs(kmObs)
!                  VENTO  - dirObs uObs(kmObs)
!                           velObs vObs(kmObs)
!
!     Para assimilar, preciso criar as variaveis contendo as coordenadas dos
!     pontos de grade. Utilizarei, tambem, as valores da variavel poderosa
!     slp(i,j,2) - pressao ao nivel medio do mar - Esquema de Messinger
!     slp(i,j,3) - pressao a superficie
!     slp(i,j,4) - altitude da superficie
!     Estes campos sao utilizados para o ajuste do geopotencial a temperatura.
!     A segunda variavel poderosa, utilizada, e datsav(1) - pontos na direcao x
!                                                     (2) - pontos na direcao y
!                                                     (3) - latitude inicial da grade
!                                                     (4) - longitude inicial da grade
!                                                     (5) - latitude final
!                                                     (6) - longitude final
!     Com as informacoes da variavel poderosa 2 vou criar as seguintes variaveis:
!     latGr(jGr),lonGr(iGr)
       latGr( 1)=float(datsav(3))/1000.
       latGr(ny)=float(datsav(5))/1000.
       deltay=(latGr(ny)-latGr(1))/float(ny-1)
       do j=2,ny
          latGr(j)=latGr(1)+float(j-1)*deltay
       enddo
       lonGr( 1)=float(datsav(4))/1000.
       lonGr(nx)=float(datsav(6))/1000.
       deltax=(lonGr(nx)-lonGr(1))/float(nx-1)
       do i=2,nx
          lonGr(i)=lonGr(1)+float(i-1)*deltax
       enddo
!       print*
!       print*,'Longitude'
!       print*,lonGr
!       print*
!       print*,'Latitude'
!       print*,latGr
!      quero identificar o nivel p para cada ponto
!      de grade
       do i=1,nx
          do j=1,ny
          k=1
          do
	       dif=slp(i,j,3)/100.0-pr(k)
	       if(dif.gt.0.) then
	          ind(i,j)=k
		  !if(ind(i,j).gt.nz) then
		  !   print*,ind(i,j),slp(i,j,3),i,j,k
		  !endif
		  exit
	       else
	          k=k+1
	       endif
	     enddo
	  enddo
       enddo             
       dadosAssim='/home2/jgomes/worketa/data/metar/assim.txt'
       ind=index(dadosAssim,' ')-1
!       dadosAssim(45:53)='assim.txt'
       open(10,file=dadosAssim,status='old',iostat=istatAss1)
       print*,dadosAssim,istatAss1
       k=1
       if(istatAss1.ne.0) then
         print*,'Problemas no arquivo a ser assimilado...'
         print*,'Nao havera assimilacao...'
       else
         NPI=4.*atan(1.0)
         do
           read(10,'(1x,3(i2,1x),a4,1x,f5.2,1x,f7.2,1x,f7.2,1x,f6.1)'
     &     ,iostat=istatAss2)
     &     dia,hora,minuto,nomeAer(k),
     &     tempObs(k),pnmmObs(k),veloc,dir
           if(istatAss2.ne.0) exit
!           write(*,'(1x,3(i2,1x),a4,4(1x,f8.2))',iostat=istatAss2)
!     &     dia,hora,minuto,nomeAer(k),
!     &     tempObs(k),pnmmObs(k),veloc,dir
           dirRad=dir*NPI/180.
           uObs(k)=-veloc*sin(dirRad)
           vObs(k)=-veloc*cos(dirRad)
	   tempObs(k)=tempObs(k)+273.16
	   pnmmObs(k)=pnmmObs(k)*100.00
           k=k+1
           if(k.gt.kmObs) then
             print*,'quantidade de pontos assimilados alcancou valor'
             print*,'acima do valor permitido kmObs ',kmObs
             exit
           endif
         enddo
         nObs=k-1
       endif
       close(10)     
       open(10,file='/home2/jgomes/worketa/data/metar/'
     &    //'locais.txt',status='old',iostat=istatAss3)
       if(istatAss3.ne.0) then
        print*,'Nao encontrei o arquivo com as coordenadas'
        print*,'das observacoes...'
        print*,'Teremos problemas...'
       else
        k=1
        do
          if(k.gt.nObs)exit
          read(10,'(a4,1x,3(f7.2,1x))',iostat=istatAss4)
     &    aero,longi,lati,alti
          if(istatAss4.ne.0) exit
          k=1
          write(*,'(a4,1x,3(f7.2,1x))',iostat=istatAss4)
     &    aero,longi,lati,alti
          do
           if(aero.eq.nomeAer(k)) then
             longObs(k)=longi
	     if(longObs(k) < 0.0) longObs(k)=360.0+longObs(k)
             latiObs(k)=lati
             altiObs(k)=alti
             !print*,nomeAer(k)
             !print*,longObs(k),latiObs(k),altiObs(k),k
             !print*,longi,lati,alti
             exit
           endif
           k=k+1
           if(k.gt.nObs)exit
          enddo
        enddo
       endif
       close(10)
       print*,'chequei aqui... Niv 1'
       print*,'entrando na ajuste dados...',nObs
       call AjusteDados(nObs,altiObs,pnmmObs,tempObs,uObs,vObs,
     &                                   slp,     tp,  uw,  vw,
     & nx,ny,nz,lonGr,latGr,longObs,latiObs,pr,ind)
!Niv
       endif
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: ',k,pr(k)
	if (datsav(11) .ne. 104) then
	     print *,'Calling RH patch.'
	     call rh_fix(nx,ny,nz,rh,pr)
	else

Cmp     moisture every 50 for grid 104...handle in special way

        if (mod(k,2) .eq. 0 .or. k.eq.nz) then

        if (k.lt.(nz-1)) then
        write(6,*) 'ave value from ', pr(K+1),pr(K-1),k+1,k-1
        DO J=1,NY
        DO I=1,NX
        RH(I,J,K)=(RH(I,J,K+1)+RH(I,J,K-1))*0.5
        ENDDO
        ENDDO
        else
        write(6,*) 'copying value from ', pr(K-1),k-1
        DO J=1,NY
        DO I=1,NX
        RH(I,J,K)=RH(I,J,K-1)
        ENDDO
        ENDDO
        endif

        endif

	endif

	  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***	Change 1000 hPa heights to match lowest temps
C***
	biassum=0.

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

C	old=ht(i,j,1)
	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
C this criteria empirically derived.  Was found that where 950-975 thick >
C about 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


  633	format(7(f6.2,1x))
	enddo
	enddo

C	call smooth(ht,nx,ny,nz,2)

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

C***
C***
C***


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 (datsav(11) .ne. 104) then
              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
	    elseif (datsav(11) .eq. 104) then
	      if (k.eq.7 .or. 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,5)
                  enddo
                enddo
              elseif(k.eq.8 .or. 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,6)
                  enddo
                enddo
              endif
	    endif
	  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
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),th(nx/2,ny/2,nz),rh(nx/2,ny/2,nz)
     .       ,uw(nx/2,ny/2,nz),vw(nx/2,ny/2,nz)
      print *,'ua :',ht(nx/2,ny/2,1),th(nx/2,ny/2,1),rh(nx/2,ny/2,1)
     .       ,uw(nx/2,ny/2,1),vw(nx/2,ny/2,1)

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

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

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

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

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

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

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

 1076   continue

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


c
c *** Create output file name.
c
      n=index(outdir,' ')-1
	if (datsav(11) .eq. 221) model='ETA221'
	if (datsav(11) .eq. 212) model='ETA212'
	if (datsav(11) .eq. 104) model='ETA104'
	if (datsav(11) .eq. 255) model='ETAwrk'
      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
C
      nsfcfld=12

	if (datsav(10) .eq. 3) gproj='LC'
	if (datsav(10) .eq. 5) gproj='PS'
	if (datsav(10) .eq. 0) gproj='CE'
	write(1) nz

        n=index(gribfile,' ')-1
	write(6,*) n
	write(6,*) gribfile
	write(6,*) 'calling get_fullgds with ',
     +    gribfile(1:n)
	call get_fullgds(gribfile(1:n),datsav(1),datsav(2),kgds)
	write(6,*) 'back from get_fullgds'

C new stuff for GDS file
C
        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')

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

	close(14)
C
C end new stuff

	write(6,*) '********************************************'
	write(6,*) 'writing out the grid with the following info'
	write(6,*) 'dimensions nx,ny,nz ', nx,ny,nz
	write(6,*) '********************************************'

c
c *** Write isobaric upper air data.
c
      write(1) uw
      write(1) vw
      write(1) ht
      write(1) rh
      write(1) pr
      write(1) slp

 	write(6,*) 'writing these sfc values to the ETA file...'
        write(6,*) 'at center of input grid'
        do k=1,nsfcfld
        write(6,*) 'field, value ', k, slp(nx/2,ny/2,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,datsav
c
      real pr(nz),ht(nxny,nz),tp(nxny,nz)
     .      ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz)
     .      ,dummy,slp(nxny,12),tmp(nxny)
      real crot(nxny),srot(nxny),rmagb,rmagaft,ubef,vbef,urlat,
     .	    urlon
c
      integer kgds(200),kpds(50),len,kerr,jpds(200),jgds(200)
     .         ,lenpds,lenkgds,nwords,kpdsl
     .        ,j,k,nx,IRETO,KNUM,IRET1
c
      character*(*) etafile
      character(LEN=1):: pds(50)
      character(LEN=12):: atime
      character(LEN=2):: gproj
      logical bitmap(nxny)

      common /gdsinfo/datsav(11)
c_______________________________________________________________________________
c
      len=index(etafile//' ',' ')-1
c
c *** Check that the input dimensions match grib file dimensions. 
c
	write(6,*) 'inside read_degrib '

CCC
CCC  Should be able to derive all header info in the DEGRIB file
CCC from the GDS/PDS of the GRIB file.
CCC
CCC
c
c *** Fill time stamp (assume all records are for same time).

	call get_gds(etafile(1:len),datsav,kpds)

	write(6,*) 'KPDS vals ', kpds(10),kpds(9),kpds(8)

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

      write(atime,'(i2.2,i2.2,i2.2,i2.2,i4.4)') 
     .   kpds(8),kpds(9),kpds(10),kpds(11),kpds(14)
c
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
Cmp
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)

	write(6,*) 'first GETGB, IRET1= ', IRET1
	write(6,*) 'LAND/SEA READ!!!!! '
	write(6,*) (slp(j,1),j=nwords/2,nwords/2+5)

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

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,2),IRET1)
        slp(:,2)=slp(:,2)*100
	write(6,*) 'PMSL READ!!!!! '
	write(6,*) (slp(j,2),j=nwords/2,nwords/2+8)

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

       call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,slp(1,3),IRET1)
        slp(:,3)=slp(:,3)*100
        write(6,*) 'PSFC READ!!!!! '
        write(6,*) (slp(j,3),j=nwords/2,nwords/2+8)

	jpds(5)=132
	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)

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)=191
	jpds(6)=1
	jpds(7)=0

       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 over ',jpds(7), 'layer!!'
        write(6,*) (slp(j,5),j=nwords/2,nwords/2+8)
	endif

	jpds(7)=0
       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 over ',jpds(7), 'layer!!', IRET1
        write(6,*) (slp(j,7),j=nwords/2,nwords/2+8)
	endif
	
	jpds(7)=0
       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 over ',jpds(7), 'layer!!',IRET1
        write(6,*) (slp(j,9),j=nwords/2,nwords/2+5)
	endif

	jpds(5)=175
        jpds(6)=1
        jpds(7)=0
       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 over ',jpds(7), 'layer!!',IRET1
        write(6,*) (slp(j,11),j=nwords/2,nwords/2+5)
	endif

C 	MOISTURE

	jpds(5)=182
	jpds(7)=0

       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 over ',jpds(7), 'layer!!', IRET1
        write(6,*) (slp(j,6),j=nwords/2,nwords/2+8)
	endif

        jpds(7)=0
       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 over ',jpds(7), 'layer!!', IRET1
        write(6,*) (slp(j,8),j=nwords/2,nwords/2+8)
	endif

        jpds(7)=0
       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 over ',jpds(7), 'layer!!',IRET1
        write(6,*) (slp(j,10),j=nwords/2,nwords/2+5)
	endif

        jpds(5)=183
	jpds(7)=0
       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 over ',jpds(7), 'layer!!',IRET1
        write(6,*) (slp(j,12),j=nwords/2,nwords/2+5)
	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))
        write(6,*) 'LEVEL ', jpds(7)
	jpds(5)=52
      call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,rh(1,k),IRET1)
	if (rh(20,k).gt.1.e20) write(6,*) 'getgbd(1) rh of ',k,rh(20,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 (rh(20,k).gt.1.e20) write(6,*) 'getgbd(2) rh of ',k,rh(20,k)
!	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+55,k),tp(nxny/2+55,k),
     +	rh(nxny/2+55,k),uw(nxny/2+55,k),vw(nxny/2+55,k)

	enddo


c
c *** Normal finish.
c
1000  continue


C	rotate the winds at this point
	nx=KGDS(2)

	write(6,*) 'using datsav(10)= ', datsav(10)
	if (datsav(10) .eq. 3) gproj='LC'
	if (datsav(10) .eq. 5) gproj='PS'

	if (gproj .eq. 'LC') then
	write(6,*) 'rotating a lambert projection'
	call rotate_lcc(kgds,crot,srot)
	elseif (gproj .eq. 'PS') then
	call rotate_str(kgds,crot,srot,urlat,urlon)
	else
	write(6,*) 'not doing a wind rotation...'
	goto 1075
	endif

	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

 1075	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
	write(6,*) 'checking RH at lev= ', pr(k)
	    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 rotate_lcc(kgds,crot,srot)

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


	parameter (NPTS=100000)
      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

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       subroutine get_gds(etafile,gdsinfo,kpds)

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

      integer kgds(200),kpds(200),len,kerr
     .         ,lenpds,lenkgds,nwords,kpdsl
     .         ,j,k,gdsinfo(11)
     .         ,gdsav,IRETO,JGDS(200),JPDS(200)
      real tmp(100000)
      logical bitmap(100000)


	nxny=100000


	JPDS=-1
	JGDS=-1


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

	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)=11
	jpds(6)=100
	jpds(7)=500
	call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS,
     &     BITMAP,tmp,IRET1)

	write(6,*) 'back from getgb in get_gds: ', 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)
	gdsinfo(10)=KGDS(1)
	gdsinfo(11)=KPDS(3)
       write(6,*) gdsinfo

	write(6,*) 'KPDSinfo ', (kpds(I),i=1,10)

        return
        print *,'GETGDS PROBLEM'
        stop

        end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

                subroutine lambert(gdslatur,gdslonur)

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

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

        common /gdsinfo/ datsav(11)

        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)

        write(6,*) 'values in lambert '
        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

        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)

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

	write(6,*) 'xll, yll, xll**2., yll**2. ', xll, yll, xll**2., 
     +	yll**2.
	write(6,*) 'lat pieces: ', rpi/2., (xll**2.+yll**2.), earth_rad,
     +	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


        gdslatur= ( rpi/2. - 2 *
     + atan ( ( (abs(xur)**2.+abs(yur)**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

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

        return

        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        subroutine rotate_str(kgds,crot,srot,urlat,urlon)

C
C Program adapted 26 March 99 for workstation Eta.  Original code from
C iplib program GDSWZD05
C
C

C SUBPROGRAM:  GDSWZD05   GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL
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 POLAR STEREOGRAPHIC AZIMUTHAL 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           OPTIONALLY, THE VECTOR ROTATIONS AND THE MAP JACOBIANS
C           FOR THIS GRID MAY BE RETURNED AS WELL.

C       POLAR STEREOGRAPHIC GRIDS
C          (2)   - N(I) NR POINTS ALONG LAT CIRCLE
C          (3)   - N(J) NR POINTS ALONG LON CIRCLE
C          (4)   - LA(1) LATITUDE OF ORIGIN
C          (5)   - LO(1) LONGITUDE OF ORIGIN
C          (6)   - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
C          (7)   - LOV GRID ORIENTATION
C          (8)   - DX - X DIRECTION INCREMENT
C          (9)   - DY - Y DIRECTION INCREMENT
C          (10)  - PROJECTION CENTER FLAG
C          (11)  - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28)

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

        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)
        H=(-1.)**IPROJ
        HI=(-1.)**ISCAN
        HJ=(-1.)**(1-JSCAN)
        DXS=DX*HI
        DYS=DY*HJ

        DE=(1.+SIN(60./DPR))*RERTH
        DR=DE*COS(RLAT1/DPR)/(1+H*SIN(RLAT1/DPR))
        XP=1-H*SIN((RLON1-ORIENT)/DPR)*DR/DXS
        YP=1+COS((RLON1-ORIENT)/DPR)*DR/DYS

        DE2=DE**2
        XMIN=0
        XMAX=IM+1
        YMIN=0
        YMAX=JM+1
        NRET=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSLATE GRID COORDINATES TO EARTH COORDINATES

        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*DPR*ATAN2(DI,-DJ)+3600,360.)
                RLAT(N)=H*DPR*ASIN((DE2-DR2)/(DE2+DR2))
              ENDIF

              NRET=NRET+1

              IF(LROT.EQ.1) THEN
                IF(IROT.EQ.1) THEN
                 CROT(N)=H*COS((RLON(N)-ORIENT)/DPR)
                 SROT(N)=SIN((RLON(N)-ORIENT)/DPR)
                ELSE
                 CROT(N)=1
                 SROT(N)=0
                ENDIF
              ENDIF
	ENDIF

          ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	urlat=rlat(im*jm)
	urlon=rlon(im*jm)
	if (urlon .gt. 180) urlon=urlon-360.
	if (urlon .lt. -180) urlon=urlon+360.
	write(6,*) 'in stereo sub found lat,lon ', urlat,urlon
C	write(6,*) 'lat/lon in stereo '
	do N=1,IM*JM
        J=INT((N-1)/IM)+1
        I=N-(J-1)*IM
C	write(6,69) I,J,RLAT(N),RLON(N)
	enddo
   69   format(I3,1x,I3,1x,2(f7.3,1x))
        RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	        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 '
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

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	subroutine smooth(a,idim,jdim,ldim,npass)
	integer idim,jdim,ldim,npass
	real a(idim,jdim,ldim)

	DO N=1,NPASS


	do L=1,LDIM
	 do J=2,JDIM-1
	  do I=2,IDIM-1

	  a(i,j,l)=1./8.*(4.*a(i,j,l)+a(i+1,j,l)+a(i-1,j,l)
     +				     +a(i,j+1,l)+a(i,j-1,l))
	  enddo
         enddo
	enddo

	ENDDO

	return
	end
!Niv
!
	subroutine AjusteDados(nObs,altO,pnmmObs,tempObs,uObs,vObs,
     &                                       slp,     tp,  uw,  vw,
     &  im,jm,km,    x,    y,    xo,      yo,pr,ind)
	implicit none
	logical contAlt
	integer i,i1,ifator,im,j,j1,jfator,jm,k,karq,km
        integer kmObs,nObs
        parameter(kmObs=100)
	integer irec
	real Rdef,tan,yVal
	real x0,y0
        integer ind(im,jm)
	real altP(im,jm),dist(im,jm,kmObs)
	real altO(kmObs)
	real W(im,jm,kmObs)
	real aux(im,jm),aux1(im,jm),aux2(im,jm)
	real x(im),xo(kmObs),y(jm),yo(kmObs)
	real pnmmObs(kmObs),tempObs(kmObs),uObs(kmObs),vObs(kmObs)
	real tp(im,jm,km),uw(im,jm,km),vw(im,jm,km)
	real pr(km)
	real slp(im,jm,12)
	do i=1,im
	   do j=1,jm
	      altP(i,j)=slp(i,j,4)
	   enddo
	enddo
	print*,'altitudes das observacoes...'
	!do k=1,nObs
	!   write(*,*) altO(k)
	!enddo
	W=9999.9
        Rdef=2.5 
        Rdef=5.0 
	Rdef=7.5
	contAlt=.true.
	if(contAlt) then ! Testando o efeito de altitude...
	do k=1,nObs
	   do j=1,jm
	     do i=1,im
	      dist(i,j,k)=sqrt((yo(k)-y(j))**2+(xo(k)-x(i))**2)
	        !write(99,*)dist(i,j,k),xo(k),yo(k),i,j,k
		!write(99,*)x(i),y(j)
		!write(99,*)altP(i,j)
		!write(99,*)
	      if(dist(i,j,k) < Rdef) then
	        if(xo(k) /= x(i)) then
	           tan=(yo(k)-y(j))/(xo(k)-x(i))
                else
                   tan=9999.99
                endif
                if((tan > 0.).and.(tan < 9999.99))then
		  if(xo(k) > x(i)) then
                     i1=i+1
		     if(i1 > im) i1=im
		     ifator=1
		     jfator=1
		  else
		     i1=i-1
		     if(i1 < 1) i1=1
		     ifator=-1
		     jfator=-1
		  endif
                  j1=j
                  do
                    !if(x(i1) < xo(k)) then
		    if((x(i1)-xo(k))*ifator < 0) then
                      yval=(x(i1)-x(i))*tan + y(j)
                      if((yval > y(j1)).and.(yval <= y(j1+1)))then
                        if((altP(i,j) < altP(i1,j1 )).and.
     &                     (altP(i,j) < altP(i1,j1+1*jfator)).and.
     &                     (altO(k) < altP(i1,j1 )).and.
     &                     (altO(k) < altP(i1,j1+1*jfator))) then
                           !O ponto i,j nao e influenciado pela obs.
                           W(i,j,k)=0.0
                           exit
			else
			   i1=i1+1*ifator
                        endif
		      else
		         j1=j1+1*jfator
		      endif
		      !if(i1 >= im) exit
		    else
		       !
		       !O ponto e influenciado.
		       !Preciso
		       W(i,j,k)=1.0
		       exit
		    endif
		  enddo
		else
		 if(tan < 0.0) then
		 !case(:0.0)
		  if(x(i) < xo(k)) then
		    i1=i+1
		    if(i1 > im) i1 = im
		    ifator=1
		    jfator=-1
		  else
                    i1=i-1
		    if(i1 < 1) i1 = 1
		    ifator=-1
		    jfator=1
		  endif
                  j1=j
		  do
		     !read*
		     !if(x(i1) > xo(k))then
		     if((x(i1)-xo(k))*ifator < 0) then
		        yval=(x(i1)-x(i))*tan+y(j)
			if((yval > y(j1)).and.(yval<=y(j1+1))) then
                        if((altP(i,j) < altP(i1,j1 )).and.
     &                     (altP(i,j) < altP(i1,j1+1*jfator)).and.
     &                     (altO(k) < altP(i1,j1 )).and.
     &                     (altO(k) < altP(i1,j1+1*jfator))) then
			     !Este ponto nao e influenciado
                             W(i,j,k)=0.0
                             exit
                          else
			    i1=i1+1*ifator
			  endif
			else
			  j1=j1+1*jfator
			  !if(j1 >= jm) then
			  !   W(i,j,k)=1.0
			  !   exit
			  !endif
			endif
			!if(i1 <=1 ) then
			!   W(i,j,k)=1.0
			!   exit
			!endif
		     else
			!O ponto e influenciado, vou calcular
			!o peso.
			W(i,j,k)=1.0
			exit
		     endif
                  enddo
		 else
		   if(tan == 0.) then
		      if(xo(k) < x(i)) then
		        j1=j
			i1=i-1
			if(i1 < 1) i1 =1
			ifator=-1
                        do
			  if((altP(i,j) < altP(i1,j1)).and.
     &                       (altO(k)   < altP(i1,j-1))) then
			     !Este ponto nao e influenciado
                             W(i,j,k)=0.0
                             exit
                          else
			    i1=i1+1*ifator
			    if((i1 < 1).or.(xo(k) > x(i1))) then
			      W(i,j,k)=1.0
			      exit
			    endif
			  endif
		        enddo
                      else
		        j1=j
			i1=i+1
			if(i1 > im) i1=im
			ifator=1
                        do
			  if((altP(i,j) < altP(i1,j1)).and.
     &                       (altO(k)   < altP(11,j1))) then
			     !Este ponto nao e influenciado
                             W(i,j,k)=0.0
                             exit
                          else
			    i1=i1+1*ifator
			    if((i1 > im).or.(xo(k) < x(i1))) then
			      W(i,j,k)=1.0
			      exit
			    endif
			  endif
		        enddo
		      endif
		   endif
		   if(tan == 9999.99) then
		      if(yo(k) < y(j)) then
		        i1=i
			j1=j-1
			if(j1 < 1) j1 =1
			jfator=-1
                        do
			  if(altP(i,j) < altP(i1,j1)) then  !.and.
			     !Este ponto nao e influenciado
                             W(i,j,k)=0.0
                             exit
                          else
			    j1=j1+1*jfator
			    if((j1 < 1).or.(yo(k) > y(j1))) then
			      W(i,j,k)=1.0
			      exit
			    endif
			  endif
		        enddo
                      else
		       if((yo(k)==y(j)).and.(xo(k)==x(i))) then
                        W(i,j,k)=1.0
		        !exit
		       else
		        i1=i
			j1=j+1
			if(j1 > jm) j1=jm
			jfator=1
                        do
			  if(altP(i,j) < altP(i1,j1)) then  !.and.
			     !Este ponto nao e influenciado
                             W(i,j,k)=0.0
                             exit
                          else
			    j1=j1+1*jfator
			    if((j1 > jm).or.(yo(k) < y(j1))) then
			      W(i,j,k)=1.0
			      exit
			    endif
			  endif
		        enddo
                       endif
		      endif
		   endif
		 endif
                endif
	      else
	        W(i,j,k)=0.0
	      endif
	     enddo
	   enddo
	enddo
	else
	   W=1.0
	endif
	print*,'vou entrar na Assimila...'
	do i=1,im
	   do j=1,jm
	      aux(i,j)=uw(i,j,ind(i,j))
              aux2(i,j)=aux(i,j)
	   enddo
	enddo
	irec=0
        call Assimila(uObs,aux,aux1,W,altP,altO,im,jm,nObs,
     &                x,y,xo,yo,irec)
	do i=1,im
	   do j=1,jm
              write(88,*) aux2(i,j)-aux1(i,j)
	      uw(i,j,ind(i,j))=aux1(i,j)
	   enddo
	enddo
	print*,'assimilei u'
        do i=1,im
           do j=1,jm
              aux(i,j)=vw(i,j,ind(i,j))
           enddo
        enddo
        call Assimila(vObs,aux,aux1,W,altP,altO,im,jm,nObs,
     &                x,y,xo,yo,irec)
        do i=1,im
           do j=1,jm
              vw(i,j,ind(i,j))=aux1(i,j)
           enddo
        enddo
	print*,'assimilei v'
	do i=1,im
	   do j=1,jm
	      aux(i,j)=tp(i,j,ind(i,j))
	   enddo
	enddo
        call Assimila(tempObs,aux,aux1,W,altP,altO,im,jm,nObs,
     &                x,y,xo,yo,irec)
	do i=1,im
	   do j=1,jm
	      tp(i,j,ind(i,j))=aux1(i,j)
	   enddo
	enddo
	print*,'assimilei Temp'
	do i=1,im
	   do j=1,jm
	      aux(i,j)=slp(i,j,2)
	   enddo
	enddo
        call Assimila(pnmmObs,aux,aux1,W,altP,altO,im,jm,nObs,
     &                x,y,xo,yo,irec)
	do i=1,im
	   do j=1,jm
	      slp(i,j,2)=aux1(i,j)
	   enddo
	enddo
	print*,'assimilei pnmm'
	print*,'saindo da Ajuste...'
        return
	end
!Niv
       subroutine Assimila(tsmObs,tsm,tsmNova,W1,altP,altO,nx,ny,nObs,
     &                     lon,lat,lonObs,latObs,irec)
        !definicao das variaveis
	implicit none
	integer irec
	integer i,istat1,istat2,istat3,im,j,jm,k,km,kmax
	integer nx,ny
        integer kmObs,nObs
        parameter(kmObs=100)
	real alt,angZen
	real delta,difLon,difLat,dist,distV
	real fatorCor,fatorCorV
	real latMax,latMin,lonMax,lonMin
	real rdef,rVert,somaPeso,somaPesoV
	real tsm1,tsm2,tsm3,tsm4
	real lon(nx),lat(ny),latObs(kmObs),lonObs(kmObs)
	real tsm(nx,ny),tsmNova(nx,ny)
	real tsmInt(kmObs),tsmObs(kmObs)
	real altP(nx,ny),altO(kmObs)
	real w1(nx,ny,kmObs)
	real w(kmObs),wVert(kmObs)
	rdef=2.5
	rdef=5.0
	rdef=7.5
	rVert=25.
	rVert=50.
	rVert=75.
	tsmNova=tsm
	latMax=-999.9
	latMin= 999.9
	lonMax=-999.9
	lonMin= 999.9
	do k=1,nObs
	     if(lonObs(k) < 0.0) lonObs(k)=360.0+lonObs(k)
	     if(latMax < latObs(k)) latMax=latObs(k)
	     if(latMin > latObs(k)) latMin=latObs(k)
	     if(lonMax < lonObs(k)) lonMax=lonObs(k)
	     if(lonMin > lonObs(k)) lonMin=lonObs(k)
	enddo
	close(10)
        !Os dados ja foram lidos...
	!Geracao das latitudes e longitudes dos campos do NMC
        do k=1,nObs
	     i=1
	     do !1
               difLon=lonObs(k)-lon(i)
	       j=1
	       if(difLon < delta) then
 	         do !2
	          difLat=latObs(k)-lat(j)
		  if(difLat < delta) then
		    tsm1=tsm(i,j)+(tsm(i+1,j)-tsm(i,j))*
     &              (lonObs(k)-lon(i))/(lon(i+1)-lon(i))
		    tsm2=tsm(i,j+1)+(tsm(i+1,j+1)-tsm(i,j+1))*
     &              (lonObs(k)-lon(i))/(lon(i+1)-lon(i))
		    tsm3=tsm(i,j)+(tsm(i,j+1)-tsm(i,j))*
     &              (latObs(k)-lat(j))/(lat(j+1)-lat(j))
		    tsm4=tsm(i+1,j)+(tsm(i+1,j+1)-tsm(i+1,j))*
     &              (latObs(k)-lat(j))/(lat(j+1)-lat(j))
	            tsmInt(k)=(tsm1+tsm2+tsm3+tsm4)*0.25
		    i=nx+1
		    exit
		  else
		    j=j+1
		  endif
		  if(j>ny)exit
	         enddo!2
	       else
	         i=i+1
	       endif
	       if(i>nx)exit
	     enddo!1
	enddo
	!read*
	!Calculo do peso das observacoe e ajustes...	
	tsmNova=tsm
        do i=1,nx
	   if((lon(i) <= lonMax).and.(lon(i)>=lonMin)) then
	      do j=1,ny
	         if((lat(j) <= latMax).and.(lat(j)>=latMin)) then
		    somaPeso=0.0
		    do k=1,nObs
		     if(W1(i,j,k)==1.0) then
		       dist=(lonObs(k)-lon(i))**2+(latObs(k)-lat(j))**2
		       dist=sqrt(dist)
		       if(dist < Rdef) then
                         w(k)=(Rdef**2-dist**2)/(Rdef**2+dist**2)
			 w(k)=exp(-1.0*dist*dist/(2.0*Rdef*Rdef))
		       else
		          w(k)=0.0
		       endif
		       distV=abs(altP(i,j)-altO(k))
!		       if(distV < RVert) then
!                        wVert(k)=(RVert**2-distV**2)/(RVert**2+distV**2)
!		       else
!		        wVert(k)=0.0
!		       endif
                       if(distV > RVert) then
                         wVert(k)=0.0
		       else
		         if(distV <= RVert/4.0) then
                           wVert(k)=1.0
			 else
			   wVert=1.0-distV/(altP(i,j)+altO(k))
			 endif
		       endif
		       somaPeso=somaPeso+w(k)
		       somaPesoV=somaPesoV+wVert(k)
		     else
		       wVert(k)=0.0
		       W(k)=0.0
		     endif
		    enddo
		    fatorCor=0.0
		    fatorCorV=0.0
!		    wVert=1.0      !SITUACAO ESPECIAL... ESTOU FAZENDO TESTES...
!		    somaPesoV=0.0  !Estou fazendo testes...
		    do k=1,nObs
		      fatorCor=fatorCor+w(k)**2*wVert(k)**2
     &                *(tsmObs(k)-tsmInt(k))
		    enddo
		    if((somaPeso+somaPesoV) > 0.0) then
		       fatorCor=fatorCor/(somaPeso+somaPesoV)
		    endif
		    tsmNova(i,j)=tsm(i,j)+fatorCor
		    write(99,'(4(e18.8,1x),2(i4,1x))')
     &              tsmNova(i,j),tsm(i,j),
     &              fatorCor,somaPeso,i,j
		 endif
	      enddo
	   endif
	enddo
	close(10)
	OPEN(10,file='dadAjus.dat',form='unformatted',status='unknown',
     &  access='direct',recl=nx*ny*4)
         irec=irec+1
         write(10,REC=irec)tsm
	 irec=irec+1
         write(10,REC=irec)tsmNova
	 irec=irec+1
         write(10,REC=irec)altP
        close(10)
        open(10,file='dados1.dat',status='unknown')
         do k=1,nObs
           write(10,'(i4,1x,2(e12.4,1x))')k,tsmObs(k),tsmInt(k)
         enddo
        close(10)
	print*,'saindo da Assimila...'
       return
       end
!Niv
