      program degrib2model_driver
c
      implicit none
c
      integer nx,ny,nz,kpds(200),datsav(9),n,k         !Degrib grid dimensions
c
      character(LEN=255):: gribfile,outfile
c
      real esat,es
c
      common /estab/esat(15000:45000),es(15000:45000)

c_______________________________________________________________________________
c
c *** Initialize tables.
c
      call es_ini
c
c *** Read file names.
c
      read(5,'(a)') gribfile
      read(5,'(a)') outfile

c!       gribfile='../data_in/init/gfsdata'//c_k//'.grb'
c!       outfile ='../data_in/init/gfsdata'//c_k//'.dat'
c
       n=index(gribfile,' ')-1
       print *,n
       write(6,*) 'first call of GET_GDS'
c      CALL GET_GDS(gribfile(1:n),datsav,kpds)!!!!!!!!!!!!!!!!!!!  DRAGAN
c      nx=datsav(1)  !!!!!!!!!!!!!!!
c      ny=datsav(2)   !!!!!!!!!!!!!!!!!!



cGSM      nx=360  !!!!!!!!!!!!!!!!!!
cGSM      ny=181  !!!!!!!!!!!!!!!!!!
      nx=1440  !!!!!!!!!!!!!!!!!!
      ny=721  !!!!!!!!!!!!!!!!!!


      if (nx .eq.  360 .and. ny .eq. 181) then
         nz=31
      elseif (nx .eq. 144 .and. ny .eq. 73) then
         nz=17
      elseif (nx .eq. 1440 .and. ny .eq. 721) then
         nz=31
      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
       print *,"nx,ny,nz,",nx,ny,nz,gribfile,outfile
c
      call degrib2model(gribfile,outfile,nx,ny,nz)

c
      end
c
c===============================================================================
c
      subroutine degrib2model(gribfile,outfile,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),ireto,n1,n2,n3,LenStr
c
cGSM      real ht(nx,ny,nz)              !Isobaric heights (m)
cGSM     .      ,tp(nx,ny,nz)              !Isobaric temps (K)
cGSM     .      ,th(nx,ny,nz)              !Isobaric theta (K)
cGSM     .      ,uw(nx,ny,nz)              !Isobaric u-wind (m/s)
cGSM     .      ,vw(nx,ny,nz)              !Isobaric v-wind (m/s)
cGSM     .      ,rh(nx,ny,nz)              !Isobaric rh,mr (%,kg/kg)
cGSM     .      ,pr(nz),pri(nz)            !Isobaric pressures (mb)
      real      lat1,lat2,lon0,sw(2),ne(2)
     .      ,xe,mrsat,esat,es
     .      ,avnlevs(31),reanlevs(17)
cGSM     .    ,slp(nx,ny,14)               !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,    dimension(:,:,:) ,allocatable    :: ht      !Isobaric heights (m)
      real,    dimension(:,:,:) ,allocatable    :: tp      !Isobaric temps (K)
      real,    dimension(:,:,:) ,allocatable    :: th      !Isobaric theta (K)
      real,    dimension(:,:,:) ,allocatable    :: uw      !Isobaric u-wind (m/s)
      real,    dimension(:,:,:) ,allocatable    :: vw      !Isobaric v-wind (m/s)
      real,    dimension(:,:,:) ,allocatable    :: rh      !Isobaric rh,mr (%,kg/kg)
      real,    dimension(:)     ,allocatable    :: pr,pri  !Isobaric pressures (mb)
      real,    dimension(:,:,:) ,allocatable    :: slp


	REAL, ALLOCATABLE:: htout(:,:,:),
     +                      uwout(:,:,:),vwout(:,:,:),rhout(:,:,:),
     +                      slpout(:,:,:)
c

      character*255 gribfile,outdir,outfile,gdsfile
      character*2   gproj
      character*11   atime
      character*8   atime_r
      character*7   model
      character*5   timex

       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, 7, 5, 3, 2, 1/

       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)


      allocate(ht(nx,ny,nz))
      allocate(tp(nx,ny,nz))                       
      allocate(th(nx,ny,nz))                     
      allocate(uw(nx,ny,nz))                 
      allocate(vw(nx,ny,nz))                
      allocate(rh(nx,ny,nz)) 
      allocate(slp(nx,ny,14)) 
      allocate(pr(nz))
      allocate(pri(nz))


c_______________________________________________________________________________
c
c *** Fill pressure levels.
c
Cmp      nsfcfld will be 12 although 4 of these are set to missing for AVN...
      print *,'degrib2model...'
      nsfcfld=14
c      nsfcfld=12    !@@@

Cmp
      write(6,*) 'using hardwired stuff! ', nx,ny,nz
         do k=1,nz
            if (nz .eq. 31) then
                pr(k)=AVNLEVS(k)
            elseif (nz .eq. 17) then
                pr(k)=reanlevs(k)
            endif
         enddo

      write(6,*) 'return read_degrib'


cGSM        OPEN(11,file='../INIT_DATA/feb1996/1996020900.bin',
cGSM     & form='unformatted', 
cGSM     & access='direct',status='old',recl=360*181*4)

        n=index(gribfile,' ')-1

	open(11,file=gribfile(1:n),form='unformatted',status='old',
     & access='direct',recl=1440*721*4)


       write(0,*) 'OPEN FILE ',gribfile(1:n)
       
c#############################################
!GSM NOVA LEITURA PARA OS ARQUIVOS GFS   

      do l=1,nz
       read(11) uw(:,:,l)
       print*,'uw l ',l,uw(241,521,l) 
      enddo

      do l=1,nz
       read(11) vw(:,:,l)
       print*,'vw l ',l, vw(241,521,l)
      enddo

      do l=1,nz
       read(11) ht(:,:,l)
       print*,' ht l ',l,ht(241,521,l)
      enddo

      do l=1,nz
       read(11) tp(:,:,l)
       print*,'tp l ',l, tp(241,521,l)
      enddo

      do l=1,nz
       read(11) rh(:,:,l)
       print*,'rh l ',l, rh(241,521,l)
      enddo

!   Missing RH at 20hPa 
!   
!GSM      rh(:,:,nz)=-99999.
!GSM      rh(:,:,nz-1)=-99999.
       

      do l=1,12
         read(11) slp(:,:,l)
         print*,'slp l ',l, slp(241,521,l)
      enddo

      close(11)


c#############################################
ccc       read(11,rec=143) ((slp(i,j,13),i=1,nx),j=ny,1,-1)
ccc       read(11,rec=144) ((slp(i,j,14),i=1,nx),j=ny,1,-1)

       
!GSM       OPEN(77,file='../INIT_DATA/feb1998/ICEC_WEASD.bin', 
       OPEN(77,file='ICEC_WEASD.bin',
     & form='unformatted',access='direct',status='old',recl=1440*721*4)
   
   
!        read(77,rec=1) ((slp(i,j,13),i=1,1440),j=721,1,-1)
!JBF teste        read(77,rec=2) ((slp(i,j,14),i=1,1440),j=721,1,-1)
    
        read(77,rec=1) ((slp(i,j,14),i=1,1440),j=1,721)
        read(77,rec=2) ((slp(i,j,13),i=1,1440),j=1,721)

      do i=1,1440
      do j=1,721
      if (slp(i,j,13).eq.9.9990003E+20) then
      slp(i,j,13)=0.
      endif
    
    
      if (slp(i,j,14).eq.9.9990003E+20) then
      slp(i,j,14)=0.
      endif
      enddo
      enddo 

        do i=1,1440
        do j=1,721

!	print*,slp(i,j,13)
        enddo
        enddo
 
!TESTE JBF:
       OPEN(888,FILE="SLP_13.txt")
       WRITE(888,*) slp(:,:,13)
       CLOSE(888)

       OPEN(889,FILE="SLP_14.txt")
       WRITE(889,*) slp(:,:,14)
       CLOSE(889)


c#########################################  
       print *,'SLP',1,slp(234,70,1)
       print *,'SLP',2,slp(234,70,2)
       print *,'SLP',3,slp(234,70,3)
       print *,'SLP',4,slp(234,70,4)

       print *,'SLP',5,slp(234,70,5)

       print *,'SLP',6,slp(234,70,6)
       print *,'SLP',7,slp(234,70,7)
       print *,'SLP',8,slp(234,70,8)
       print *,'SLP',9,slp(234,70,9)
       print *,'SLP',10,slp(234,70,10)
       print *,'SLP',11,slp(234,70,11)
       print *,'SLP',12,slp(234,70,12)

       print *,'SLP',13,slp(234,70,13)
       print *,'SLP',14,slp(234,70,14)
   
       do l=5,12
c!       print *,'SLP',14,slp(:,:,l)
       enddo
   
c       print*,slp(:,:,5)
       
       
c       print*,rh(22,:,23)
       
       
       close(11)
      close(77)
ccc      close(12)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      print *,'rh(1,1,1),',rh(1,1,1)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!GSM
!GSM      do i=1,nx
!GSM      do j=1,ny
  
!GSM        if (slp(i,j,1)<0.1) then
!GSM          slp(i,j,1)=0
!GSM        else
!GSM          slp(i,j,1)=1    
!GSM        endif        
!GSM      enddo
!GSM      enddo 


      do i=1,nx
      do j=ny,1,-1

      do l=1,nz
      if (ht(i,j,l).eq.9.9990003E+20) ht(i,j,l)=-99999.
      if (uw(i,j,l).eq.9.9990003E+20) uw(i,j,l)=-99999.
      if (vw(i,j,l).eq.9.9990003E+20) vw(i,j,l)=-99999.
      if (tp(i,j,l).eq.9.9990003E+20) tp(i,j,l)=-99999.
      if (rh(i,j,l).eq.9.9990003E+20) rh(i,j,l)=-99999.
      enddo

!GSM
      do l=1,12
      if (slp(i,j,l).eq.9.9990003E+20) slp(i,j,l)=0.
      enddo
!GSM

c############################################## DRAGAN      
      
      do l=1,21
          if (rh(i,j,l).gt.100)  rh(i,j,l)=100
      enddo
      
      do l=1,nz
      if (rh(i,j,l).lt.0) rh(i,j,l)=1
      enddo
      
      enddo
      enddo
c###########################################################################      

       print*,slp(360,:,3) 

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: ',k,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),uw(1,1,k)
            stop
         elseif (vw(1,1,k) .eq. -99999.) then
            print *,'V-wind data missing at level: ',pr(k)
            stop
         endif
      enddo

      do k=1,nsfcfld
           if (slp(1,1,k) .eq. 9.9990003E+20) then
              print *,'SLP data missing at level: ',pr(k)
           endif
      enddo
Cmp
Cmp     Handle missing surface data....
Cmp

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

CIf 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.and.k.ne.2) then  !!!!!!!!!!!?????????????????
               slp(i,j,k)=0.14
            else
               slp(i,j,k)=275.
            endif
         enddo
         enddo
         enddo
      endif

!!!	do k=5,nsfcfld

        do k=5,12

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

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

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

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

          else

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

          endif
Cmp
        enddo

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

CC!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
c!	if (nz .eq. 17) goto 969
      
cgsm      do k=1,nz

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

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

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

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

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

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

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

cgsm         if (k.eq.1) write(6,*) 'flipping isobaric data - '

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

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

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

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

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


Cmp    generalize eventually!   ???????????????????????????????????????????????????
cgsm      do k=1,nsfcfld

C       write(6,*) 'flipping surface fields... ' 

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

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

cgsm      enddo
CC!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
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,14))

          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(outfile,' ')-1

Cmp  make ETA_avn to differentiate from files based on eta grib input
      if (nz .eq. 31) then
         model='ETA_avn'
      elseif (nz .eq. 17) then
         model='ETA_rnl'
         nx=145
      endif
      n=index(outfile,' ')-1
      
cccccccccccccc      outfile=outdir(1:n)//'/'//atime(1:LenStr)//'.'//model
      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...
c      nsfcfld=12   !@@@
       nsfcfld=14
 
      gproj='LL'
      write(1) nz

C new stuff for GDS file
C

C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

C new stuff for GDS file
C
        n=index(gribfile,' ')-1
        write(6,*) 'calling get_fullgds with ',
     +    gribfile(1:n)
c       call get_fullgds(gribfile(1:n),nx,ny,kgds)
cGSM        kgds(2)=360
cGSM        kgds(3)=181
        kgds(2)=1440
        kgds(3)=721

        kgds(4)=-90000
      
        kgds(6)=128
        kgds(7)=90000
       
        kgds(8)=-1000
        kgds(9)=1000
        kgds(10)=1000
       
        kgds(20)=255



c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        write(6,*) 'back from get_fullgds'
c       call baclose(11,ireto)


        n=index(outfile,' ')-1
        gdsfile='gdsinfo.ETA_avn'

        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. 31) 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(1,1,k)
        enddo
cGSM     do k=1,nsfcfld
cGSM       write(6,*) 'field, value ', k, slp(81,261,k)
cGSM     enddo

c
      close(1)
c
      return
      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
         do j=1,ny
           do i=1,nx
!	         print *,i,j,k,rh(i,j,k)
           enddo
         enddo

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

