subroutine pusi(spl,htll,qtll,utll,vtll,ldm,pd,t,q,u,v,fis,im,jm,nm,lm)

      use dynam

      implicit none

      integer, intent(in) :: im,jm,nm,lm

      include 'const.h'
!      include 'dynam_comm.h' !!!!
      
      real,parameter::g=9.80616,gor=9.80616/287.04,epsq=1.e-12
      integer::ldm
      real,dimension(ldm)::spl,zsl,rzsl,rhsl
      real,dimension(ldm+1)::y2,pp,qq,zslh,uij,vij,qij,hsp
      real,dimension(lm+1)::sg,dg,zh
      real,dimension(lm)::dsg,sgml,qtil,zqtil,util,vtil,zuv
      real,dimension(0:im+1,0:jm+1,nm)::hgt,sm,pd,fis
      real,dimension(0:im+1,0:jm+1,nm,ldm)::htll,qtll,utll,vtll
      real,dimension(0:im+1,0:jm+1,nm,lm+1)::z,alpi
      real,dimension(0:im+1,0:jm+1,nm,lm)::t,q,u,v
!!!!!!!!!!!       
!      real,dimension(lp1) :: eta
!!!!!!!!!!!      
      real::alpt,ztop,alp,zld0,zld1,zld2,hgtp,hld0,hld2,hld1,d1,d2, &
            x,zsp,plpi,pdp,zl,zu,qld2,zss,qld0,qld1,qp,zslpu,zus,zuw, &
            zue,zun,zls,zlw,zle,zln
      integer::l,ld,i,j,n,k,lold,lp1
      integer,parameter::nsmud=1
      real,dimension(0:im+1,0:jm+1,nm)::hgts
!!      lp1=lm+1
      sg=eta

      do l=1,lm
        dsg(l)=sg(l+1)-sg(l)
      enddo

      do l=1,lm
        sgml(l)=0.5*(sg(l)+sg(l+1))
!        write(6,*) 'L, sgml(L): ', L,sgml(L)
      enddo

!-----------------------------------------------------------------------
       do ld=1,ldm+1
         y2(ld)=0.
       enddo
!--------------pusi derived constants ----------------------------------
      alpt=log(pt)
      ztop=alpt*alpt

      do ld=1,ldm
        alp=log(spl(ld))
        zsl(ld)=alp*alp
      enddo
!-----------------------------------------------------------------------
      zld2=zsl(ldm-2)
      zld1=zsl(ldm-1)
      zld0=zsl(ldm  )

      open(1,file='topo.dat', status='old',form='unformatted')
      read(1) hgt,sm
      close(1)

      call bocoh(hgt,im,jm,nm)
      call bocoh(sm,im,jm,nm)
!--------------5-point smoothing of mountains---------------------------
      if(nsmud.gt.0)    then
!-----------------------------------------------------------------------
         do k=1,nsmud
!-----------------------------------------------------------------------
         do n=1,nm
         do j=1,jm
         do i=1,im
          if(sm(i,j,n).lt.0.5)    then
             hgts(i,j,n)=(hgt(i,j-1,n)+hgt(i,j+1,n)  &
               +hgt(i+1,j,n)+hgt(i-1,j,n)          &
               +hgt(i,j,n)*4.)*0.125
          else
!             hgts(i,j,n)=hgt(i,j,n)
             hgts(i,j,n)=0.
          endif
         enddo
         enddo
!
         do j=1,jm
         do i=1,im
           hgt(i,j,n)=hgts(i,j,n)
         enddo
         enddo
!-----------------------------------------------------------------------
         enddo
	 call bocoh(hgt,im,jm,nm)
         enddo
!-----------------------------------------------------------------------
       endif

!--------------computation of pd, zeta and related variables------------
      do n=1,nm
      do j=1,jm
      do i=1,im
!-----------------------------------------------------------------------
        hgtp=hgt(i,j,n)
        hld0=htll(i,j,n,ldm)
!-----------------------------------------------------------------------
        if(hld0.gt.hgtp)    then
!--------------if sfc below lowest p sfc, extrapolate downwards---------
          hld2=htll(i,j,n,ldm-2)
          hld1=htll(i,j,n,ldm-1)
!
          d1=(zld0-zld1)/(hld0-hld1)
          d2=0.
          x=hgtp-hld0
          zsp=d2*x*x+d1*x+zld0
!-----------------------------------------------------------------------
        else
!--------------otherwise, use spline to interpolate---------------------
          do ld=1,ldm
            rhsl(ld)=htll(i,j,n,ldm+1-ld)
            rzsl(ld)=zsl(ldm+1-ld)
            y2(ld)=0.
          enddo
!
          call spline_pus(ldm,rhsl,rzsl,y2,1,hgtp,zsp,pp,qq)
!-----------------------------------------------------------------------
        endif
!-----------------------------------------------------------------------
        z(i,j,n,lm+1)=zsp
        plpi=sqrt(zsp)
        alpi(i,j,n,lm+1)=plpi
        pdp=exp(plpi)-pt
        pd(i,j,n)=pdp
!
        alpi(i,j,n,1)=alpt
        z(i,j,n,1)=ztop
!
        do l=2,lm
          plpi=log(pt+pdp*sg(l))
          alpi(i,j,n,l)=plpi
          z(i,j,n,l)=plpi*plpi
        enddo
!-----------------------------------------------------------------------
       enddo
       enddo
       enddo

!--------------velocity spline interpolation inside the domain----------
       do n=1,nm
       do j=1,jm-1
       do i=1,im-1
       do ld=1,ldm
         zslh(ld)=zsl(ld)
         uij(ld)=utll(i,j,n,ld)
         vij(ld)=vtll(i,j,n,ld)
        enddo
!
         zus=ztop
         zuw=ztop
         zue=ztop
         zun=ztop
!
        do l=1,lm
          zls=z(i,j,n,l+1)
          zlw=z(i,j+1,n,l+1)
          zle=z(i+1,j,n,l+1)
          zln=z(i+1,j+1,n,l+1)
!
          zuv(l)=.125*(zus+zls+zuw+zlw+zue+zle+zun+zln)
!
          zus=zls
          zuw=zlw
          zue=zle
          zun=zln
        enddo
!
        zslpu=(z(i,j,n,lm+1)+z(i+1,j,n,lm+1)   &
           +z(i,j+1,n,lm+1)+z(i+1,j+1,n,lm+1))*0.25
!
        if(zslpu.le.(zld0+1.e-4))then
          lold=ldm
        else
          lold=ldm+1
          zslh(lold)=zslpu
          uij(lold)=uij(lold-1)
          vij(lold)=vij(lold-1)
        endif
!
        do ld=1,lold
          y2(ld)=0.
        enddo
        call spline_pus(lold,zslh,uij,y2,lm,zuv,util,pp,qq)
        call spline_pus(lold,zslh,vij,y2,lm,zuv,vtil,pp,qq)

        do l=1,lm
          u(i,j,n,l)=util(l)
          v(i,j,n,l)=vtil(l)
        enddo
      enddo
      enddo
      enddo
!--------------computation of sigma temperatures------------------------
      do n=1,nm
      do j=1,jm
      do i=1,im
!-----------------------------------------------------------------------
       zh(1)=ztop
!
       do l=2,lm+1
        zh(l)=z(i,j,n,l)
       enddo
!
       do ld=1,ldm
         zslh(ld)=zsl(ld)
         hsp(ld)=htll(i,j,n,ld)
       enddo
!
       if(z(i,j,n,lm+1).le.(zld0+1.e-4)) then
         lold=ldm
        else
         lold=ldm+1

         zslh(lold)=z(i,j,n,lm+1)
         hsp(lold)=hgt(i,j,n)
        endif

        do ld=1,lold
          y2(ld)=0.
        enddo
!--------------temperatures inside the integration domain---------------
        call spline_pus(lold,zslh,hsp,y2,lm+1,zh,dg,pp,qq)
!
        do l=1,lm
          t(i,j,n,l)=(dg(l)-dg(l+1))*gor/(alpi(i,j,n,l+1)-alpi(i,j,n,l))
        enddo
!-----------------------------------------------------------------------
        enddo
        enddo
        enddo
!-----------------------------------------------------------------------

!--------------spec hum spline interpolation inside the domain----------
         do n=1,nm
         do j=1,jm
         do i=1,im
           do ld=1,ldm
             zslh(ld)=zsl(ld)
             qij(ld)=qtll(i,j,n,ld)
           enddo
!
           zu=ztop

           do l=1,lm
             zl=z(i,j,n,l+1)
             zqtil(l)=.5*(zu+zl)
             zu=zl
           enddo
!
          if(z(i,j,n,lm+1).le.(zld0+1.e-4))    then
            lold=ldm
          else
            lold=ldm+1
!
            zss=z(i,j,n,lm+1)
            zslh(lold)=zss
!
            qld2=qij(ldm-2)
            qld1=qij(ldm-1)
            qld0=qij(ldm  )
!
            d1=(qld0-qld1)/(zld0-zld1)
            d2=0.
            x=zss-zld0
!
            qij(lold)=d2*x*x+d1*x+qld0
          endif

          do ld=1,lold
            y2(ld)=0.
          enddo
!
          call spline_pus(lold,zslh,qij,y2,lm,zqtil,qtil,pp,qq)
!
          do l=1,lm
            qp=amax1(qtil(l),epsq)
            t(i,j,n,l)=t(i,j,n,l)/(qp*0.608+1.)
            q(i,j,n,l)=qp
          enddo
          enddo
          enddo
          enddo

	  fis=hgt*g

          call bocoh(fis,im,jm,nm)

      end subroutine pusi