!===============================================================================
!
      subroutine table(ptbl,ttbl,pt  &
                      ,rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)

use ieee_arithmetic !JBF

!     implicit none
!
! *** Generate values for look-up tables used in convection.
!
      integer,parameter::itb=76,jtb=134 !JBF: D. Chou, valores iguais do Eta regional. 
      real:: thh,ph,pq0,a1,a2,a3,a4,r,cp,eliwv,eps
!mp      parameter (thh=350.,ph=105000.
!       changed 3/8/99 to make more like operation grdeta stuff
!       appeared to have no change on underflow problem
!JBF  parameter (thh=365.,ph=105000.  &
      parameter (thh=365.,ph=105000.  &
                ,pq0=379.90516  &
                ,a1=610.78,a2=17.2693882,a3=273.16,a4=35.86  &
                ,r=287.04,cp=1004.6,eliwv=2.683e6,eps=1.e-10) !JBF: D. Chou, o eps que o Jorge comentou. 
!
!      real:: ptbl(itb,jtb),ttbl(jtb,itb),qsold (jtb),pold(jtb)  & ...........told(jtb),theold(jtb),the0 (itb),sthe(itb)  &
      real:: ptbl(itb,jtb),ttbl(jtb,itb),pold(jtb)  &
            ,qs0 (jtb),sqs   (jtb),qsnew(jtb)  &
            ,y2p (jtb),app   (jtb),aqp  (jtb),pnew(jtb)  &
            ,told(jtb),the0 (itb),sthe(itb)  &
            ,y2t (jtb),thenew(jtb),apt  (jtb),aqt (jtb),tnew(jtb)

!      double precision :: sqsk, qsold(jtb), theold(jtb)
      real :: sqsk, qsold(jtb), theold(jtb)

!_______________________________________________________________________________
 open (126,file= 'table_input.txt') !JBF
 write (126,*) 'PTBL=',ptbl,'TTBL=',ttbl,'PT=',pt,'RDQ=',rdq,'RDTH=',rdth,'RDP=',rdp,'RDTHE=',rdthe,'PL=',pl,'THL=',thl,'QS0=',qs0,'SQS=',sqs,'STHE=',sthe,'THE0=',the0
 close(126)
!
! *** Coarse look-up table for saturation point.
!
      kthm=jtb
      kpm=itb
      kthm1=kthm-1
      kpm1=kpm-1
!
      pl=pt
!
      dth=(thh-thl)/real(kthm-1)
      dp =(ph -pl )/real(kpm -1)
!
      rdth=1./dth
      rdp=1./dp
      rdq=kpm-1
!
      th=thl-dth
!

!
     open (125,file='Err.txt') !JBF
     open(128,File="err_1_laco.txt")
      do kth=1,kthm
         th=th+dth
         p=pl-dp


         do kp=1,kpm
 
            p=p+dp
            ape=(100000./p)**(r/cp)
            DENOM=th-a4*ape
            if (DENOM .GT. eps) then
             qsold(kp)=pq0/p*exp(a2*(th-a3*ape)/DENOM)
            else 
             qsold(kp)=0.
            endif

!            if (ape  .gt. 3.320840) &
!               ape=3.320840

!            qsold(kp)=pq0/p*exp(a2*(th-a3*ape)/(th-a4*ape))

            pold(kp)=p

        
        write(128,*) "kp=",kp,"  p=",p,"  ape=",ape,"  th=",th,"  qsold(kp)=",qsold(kp)
   
         enddo

!
         qs0k=qsold(1)
         sqsk=qsold(kpm)-qsold(1)
         qsold(1  )=0.
         qsold(kpm)=1.
!
         do kp=2,kpm1
            open(127,File="table.txt") !JBF
           
         if (ieee_is_finite(qs0k) .or. ieee_is_finite(sqsk)) then
         else
         write(127,*) "kp=",kp,"qsold(kp)=",qsold(kp),"qs0k=",qs0k,"sqsk=",sqsk !JBF
         endif
         close(127)
            qsold(kp)=(qsold(kp)-qs0k)/sqsk
!
! ********* Fix due to cyber half prec. limitation.
!
            if (qsold(kp)-qsold(kp-1) .lt. eps)  &
               qsold(kp)=qsold(kp-1)+eps
         if (ieee_is_NaN(qsold(kp))) write (125,*) kp
         write (125,*) pold(kp),ape,th,qsold(kp) !JBF


! ********* End fix.
!
         enddo
!
         qs0(kth)=qs0k
         sqs(kth)=sqsk
         qsnew(1  )=0.
         qsnew(kpm)=1.
         dqs=1./real(kpm-1)
!
         do kp=2,kpm1
            qsnew(kp)=qsnew(kp-1)+dqs
         enddo
!
         y2p(1   )=0.
         y2p(kpm )=0.
!
         call spline(jtb,kpm,qsold,pold,y2p,kpm,qsnew,pnew,app,aqp)
!
         do kp=1,kpm
            ptbl(kp,kth)=pnew(kp)
         enddo
!
      enddo
    close (125)
    close(128)
! *** Coarse look-up table for t(p) from constant the.
!
      p=pl-dp
      do kp=1,kpm
         p=p+dp
         th=thl-dth
         do kth=1,kthm
            th=th+dth
            ape=(100000./p)**(r/cp)
!JBF: Modificacoes

            DENOM=th-a4*ape
            if (DENOM .GT. eps) then
             qs=pq0/p*exp(a2*(th-a3*ape)/DENOM)
            else 
             qs=0.
            endif
!            qs=pq0/p*exp(a2*(th-a3*ape)/(th-a4*ape))
            told(kth)=th/ape
            theold(kth)=th*exp(eliwv*qs/(cp*told(kth)))
         enddo
!
         the0k=theold(1)
         sthek=theold(kthm)-theold(1)
         theold(1   )=0.
         theold(kthm)=1.
!
         do kth=2,kthm1
            theold(kth)=(theold(kth)-the0k)/sthek
!
            if (theold(kth)-theold(kth-1) .lt. eps)  &
               theold(kth)=theold(kth-1)+eps
!
         enddo
!
         the0(kp)=the0k
         sthe(kp)=sthek
         thenew(1  )=0.
         thenew(kthm)=1.
         dthe=1./real(kthm-1)
         rdthe=1./dthe
!
         do kth=2,kthm1
            thenew(kth)=thenew(kth-1)+dthe
         enddo
!
         y2t(1   )=0.
         y2t(kthm)=0.
!
         call spline(jtb,kthm,theold,told,y2t,kthm,thenew,tnew,apt,aqt)
!
         do kth=1,kthm
            ttbl(kth,kp)=tnew(kth)
         enddo
!
      enddo
!
      return
      end

