      subroutine vinterp(prin,h,mr,u,v,nzin,pd,t,mreta,uetat,vetat,im,jm,nm,lm)

      use dynam

      implicit none

      integer, intent(in) :: im,jm,nm,lm
!!!      include 'dynam_comm.h'
      include 'const.h'
      


      integer::nzin,nzinm1,ld,l,l1,l2,l3,i,j,n,ivi,lp1
      real,dimension(nzin)::prin
      real,dimension(0:im+1,0:jm+1,nm,nzin)::h,u,v,mr
      real,dimension(0:im+1,0:jm+1,nm)::pd
      real,dimension(0:im+1,0:jm+1,nm,lm)::t,mreta,uetat,vetat
      integer,dimension(lm)::ldt1
      real,dimension(lm)::dlt
      real::prinm,rprin,h1b,h2,h3,alpgt,hetaij,alpeta  &
           ,tetamax,alpetk,uld,vld,cf,ukb,vkb,phub,phlb,qseta, &
           rheta
      real,dimension(nzin)::alpr,alp,alpsq
      real,allocatable,dimension(:,:,:,:)::heta
      real,allocatable,dimension(:,:,:)::pdvp
!!!!!!!!!      
!      real,dimension(lm) :: deta, aeta
!      real,dimension(lp1) :: eta
!!!!!!!!!      
      integer,dimension(0:im+1,0:jm+1,nm)::lhgt
      real,dimension(0:im+1,0:jm+1,nm)::hgt,ref
      real,parameter::g=9.80616,r=287.04,pq0=379.90516,a2=17.2693882, &
              a3=273.16,a4=35.86,tresh=0.95
      real::alp1,alp2,pres

!!!      lp1=lm+1


      print *,'vinterp start'

      open(1,file='forVt.dat',status='old',form='unformatted')
      read(1)lhgt,hgt,ref
      close(1)

      call etaPr(prin,eta,pt,ldt1,nzin,lm)

      print *,'Vertically interpolate data to ETA grid.'
      print *,'input levels, output levels', nzin,LM
      print *,'in interp....pt= ', pt
!
!c
!-----------------------------------------------------------------------
!c
!c *** Vertical interpolation: quadratic interpolation of
!c        heights (equivalent to temperature being linear in ln p)
!c        and linear interpolation of winds.
!c
!c *** Computation of the 'sea level' pressure difference.
!c
      nzinm1=nzin-1
      prinm=prin(nzin)
      rprin=prinm/prin(nzinm1)
      do ld=1,nzinm1
         alpr(ld)=alog(prin(ld+1)/prin(ld))
      enddo

      do ld=1,nzin
         alp(ld)=alog(prin(ld))
         alpsq(ld)=alp(ld)**2
      enddo

      do l=1,lm
         l1=ldt1(l)
         dlt(l)=alpr(l1+1)*(alp(l1+2)-alp(l1))*(-alpr(l1))
      enddo

!Cmp *************** START FIRST SET OF LOOPS ********************

       do n=1,nm
       do j=1,jm
       do i=1,im

            l=min(lm,lhgt(i,j,n))
            l1=ldt1(l)
            l2=l1+1
            l3=l1+2
            h1b=h(i,j,n,l1)
            h2 =h(i,j,n,l2)
            h3 =h(i,j,n,l3)
            call quadInterpLnPr(h1b,h2,h3,alp(l1),alp(l2),alp(l3)  &
                           ,hgt(i,j,n),alpgt)
            pd(i,j,n)=(exp(alpgt)-pt)/ref(i,j,n)
      enddo
      enddo
      enddo

        ALLOCATE(HETA(0:IM+1,0:JM+1,nm,LM))

        do l=lm,1,-1
         do n=1,nm
         do i=1,im
         do j=1,jm


            l1=ldt1(l)
            l2=l1+1
            l3=l1+2
            h1b=h(i,j,n,l1)
            h2 =h(i,j,n,l2)
            h3 =h(i,j,n,l3)

            alpeta=alog(pt+pd(i,j,n)*eta(l))
            call quadInterpH(h1b,h2,h3,alp(l1),alp(l2),alp(l3) &
                   ,hetaij,alpeta)
            heta(i,j,n,l)=hetaij

            h1b=mr(i,j,n,l1)
            h2 =mr(i,j,n,l2)
            h3 =mr(i,j,n,l3)
            alpeta=alog(pt+pd(i,j,n)*aeta(l))
            call quadInterpH(h1b,h2,h3,alp(l1),alp(l2),alp(l3) &
                   ,hetaij,alpeta)
            mreta(i,j,n,l)=hetaij

        enddo
        enddo
        enddo
        enddo

        ALLOCATE(PDVP(0:IM+1,0:JM+1,nm))

      do n=1,nm
      DO J=1,JM-1
      DO I=1,IM-1
          PDVP(I,J,n)=0.25*(PD(I,J,n)+PD(I+1,J,n)  &
                      +PD(I,J+1,n)+PD(I+1,J+1,n))
      ENDDO
      ENDDO
      enddo

!c *** Ground surface heights converted to geopotentials, and 
!c        geopotential to temperature inversion.
!c

      tetamax = 0
      do n=1,nm
      do j=1,jm
      do i=1,im

         phub=0.
         do ivi=1,lm
            l=lm+1-ivi
            phlb=phub
          phub=g*heta(i,j,n,l)
!JBF: ----------------------------------------------------
!            write(*,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l)
!!           if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta


!!!CHOU           t(i,j,n,l)=-(phlb-phub)*(pt+aeta(l)*pd(i,j,n)) &
!!!CHOU                      /(r*deta(l)*pd(i,j,n))  !!??
!JBF:
!!           if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "i=",i, "j=",j, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta
!!!!           if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "n=", n, "i=",i, "j=",j, "aeta(l)=", aeta(l), "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l) 

!!!!!!!!!!!!1
            alp1=alog(pt+pd(i,j,n)*eta(l+1))
            alp2=alog(pt+pd(i,j,n)*eta(l))
            t(i,j,n,l)=-(phlb-phub)/(r*(alp1-alp2)) !CHOU
            if ( l == lm ) t(i,j,n,l) = 200/(r*(alp1-alp2)) !CHOU

!           if (t(i,j,n,l) > 320 ) write(210,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2  

!GSM           if (i == 123 .and. j == 377 .and. n == 2)  write(210,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l),"aeta(l)=", aeta(l), &
!GSM          & "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2, "phlb-phub=", phlb-phub, "hgt(i,j,n)=", hgt(i,j,n)


            if (t(i,j,n,l).gt.tetamax) tetamax=t(i,j,n,l)

!!!!!!!!!!!!teste ANDRE E JOAO
!Cmp	substitute the TETA values at lowest level?
!	if (ivi .eq. lm .and. t(i,j,n,lm) > 320) then	
!	t(i,j,n,lm)=t(i,j,n,lm-1)
!	endif
!        if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "i=",i, "j=",j, "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l), "t(i,j,n,l)=", t(i,j,n,l)
!!!!!!!!!!!teste ANDRE E JOAO


            t(i,j,n,l)=amin1(t(i,j,n,l),325.)
            t(i,j,n,l)=amax1(t(i,j,n,l),150.)

!           if (t(i,j,n,l) > 320 ) write(220,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2  


!GSM           if (i == 123 .and. j == 377 .and. n == 2)  write(220,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l),"aeta(l)=", aeta(l), &
!GSM          & "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2, "phlb-phub=", phlb-phub, "hgt(i,j,n)=", hgt(i,j,n)  
!JBF
!            if (t(i,j,n,1) .gt. 248.3767) then  
!            t(i,j,n,1)=248.3767
!            endif

            if (eta(l+1) .le. ref(i,j,n)) then
!             pres=exp(0.5*(alp1+alp2))
              qseta=pq0/(pd(i,j,n)*aeta(l)+pt)  &
!             qseta=pq0/(pres)  &
                         *exp(a2*(t(i,j,n,l)-a3)/(t(i,j,n,l)-a4))
!JBF: ----------------------------------------------------
!            if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta
!
              rheta=mreta(i,j,n,l)/qseta
            endif
            rheta=min(rheta,tresh)
!!!!!!!!!!            mreta(i,j,n,l)=rheta*qseta  !!??
            mreta(i,j,n,l)=max(0.,mreta(i,j,n,l))
!!!!!!!!!!            t(i,j,n,l)=t(i,j,n,l)/(mreta(i,j,n,l)*0.61+1) !!??

          enddo             
!c
!C************************
!c ************ Now redefine pd to have it equal to ps-pt.
!C**************************
!c
               pd(i,j,n)=ref(i,j,n)*pd(i,j,n)
      enddo
      enddo
      enddo

!JBF
!!!     close(200)
!JBF
      DEALLOCATE(HETA)

      DO 400 L=1,LM
      DO n=1,nm
      do j=1,jm
      do i=1,im

       
        ALPETK=ALOG(PT+PDVP(I,J,n)*AETA(L))
        ld=2
        do while (ALPETK .GT. ALP(LD) .AND. LD .LT. nzin)
          ld=ld+1
        enddo
          ULD=U(I,J,n,LD)
          VLD=V(I,J,n,LD)
          CF=(ALP(LD)-ALPETK)/ALPR(LD-1)
          UKB=ULD+(U(I,J,n,LD-1)-ULD)*CF
          VKB=VLD+(V(I,J,n,LD-1)-VLD)*CF

          UETAT(I,J,n,L)=UKB
          VETAT(I,J,n,L)=VKB
        ENDDO
        enddo
        enddo
  400 CONTINUE

        DEALLOCATE(PDVP)

!c
! *** Set u, v  equal to zero at points below the ground.
!c
             do l=1,lm
             do n=1,nm
             do j=1,jm-1
             do i=1,im-1
               if (eta(l+1) .gt. ref(i,j,n)) then
                  uetat(i-1,j-1,n,l)=0.
                  vetat(i-1,j-1,n,l)=0.
                  uetat(i-1,j,n,l)=0.
                  vetat(i-1,j,n,l)=0.
                  uetat(i,j-1,n,l)=0.
                  vetat(i,j-1,n,l)=0.
                  uetat(i,j,n,l)=0.
                  vetat(i,j,n,l)=0.
               endif
            enddo
            enddo
            enddo
            enddo

      print *,"finish vinterp"
      END SUBROUTINE vinterp

