c    @(#)thermvt.f	1.1
c
c
c
c
c
c
        subroutine thermvt(psi,tvort,kz,msq)
      include'qg_par.inc'
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)
        dimension tvort(mn0,k0),psi(mn0,k0),vtm(k0,k0)
c
       do 5 k=1,kz
       do 6 kk=1,kz
       vtm(k,kk)=0.0
 6    continue
      do 7 ip=1,msq
      tvort(ip,k)=0.
 7    continue
 5    continue
c
c
      do 20 i=1,kz
        ip=0
        do 19 j=1,kz
          tmp=0.0
          do 18 k=1,kz
            tmp=tmp+amat(i,k)*eigval(k)*ainv(k,j)
   18     continue
          vtm(i,j)=tmp
   19   continue
   20 continue
c
c --- now calculate the thermal vorticity
c
       do 500 k=1,kz
        do 170 kk=1,kz
          tmp=vtm(k,kk)
          do 165 ip=1,msq
            tvort(ip,k)=tvort(ip,k)+tmp*psi(ip,kk)
  165     continue
  170   continue
 500    continue
c
c
        return
        end
