subroutine stresstkehurr(xvuse,il,iu,jl,ju,lls,nvis,iusetest)
! this subroutine determines momentum stress 
      use gridsetup
      use arrays
      use metryic
      use constants
      use restarta
      use newton
      use msga
      double precision xvuse(il:iu,jl:ju,lls,nvis)
      double precision,allocatable:: vvec(:,:,:,:), &
      rhot(:,:,:),tke(:,:,:),tkerho(:,:,:)
      double precision::kmtensor(3,3)
      double precision::velzp(3),velzm(3),velxp(3),velxm(3),velyp(3),velym(3)
      double precision,allocatable::tautilda(:,:,:,:),tauij(:,:,:,:)
      double precision::dr0(3),du0(3,3),xjac0(3,3)
      logical::made_tmp_zero
      dr0=0.5*(/ dx,dy,dz /)

      rkm=0.
      constantdiff=0.0*2.5*2.5
      if(itke.eq.1) ftke=0.
      kmtensor(1,1)=dx
      kmtensor(2,2)=dy
      kmtensor(3,3)=dz

! load model data into this subroutine
      allocate (vvec(1-ih:np+ih,1-ih:mp+ih,l,3))
      allocate (rhot(1-ih:np+ih,1-ih:mp+ih,l))
      allocate (tke(1-ih:np+ih,1-ih:mp+ih,l))
      allocate (tkerho(1-ih:np+ih,1-ih:mp+ih,l))
      allocate(tautilda(1-ih:np+ih+1,1-ih:mp+ih+1,l+1,3))
      allocate(tauij(1-ih:np+ih+1,1-ih:mp+ih+1,l+1,3))

      if(idust.eq.0) ifco2=1
      if(idust.eq.1) ifco2=2
      do 100 isweep=1,1
      ims=1
      imf=3
      if(isweep.eq.2) ims=5
      if(isweep.eq.2) imf=7
! isweep goes over both momentum variables
      do k=1,l
        do j=1-ih*j3,mp+ih*j3
          do i=1-ih,np+ih
            vvec(i,j,k,1:3)=xvuse(i,j,k,ims:imf)/xvuse(i,j,k,nv)
            rhot(i,j,k)=xvuse(i,j,k,nv)
            if(ihurrmulti.eq.1) &
            rhot(i,j,k)=xvuse(i,j,k,nv)/volgas(i,j,k)
            if(itke.eq.1)&
            tke(i,j,k)=xvuse(i,j,k,nv-1)/xvuse(i,j,k,nv)
            if(itke.eq.1)&
            tkerho(i,j,k)=xvuse(i,j,k,nv-1)
            if(itke.eq.1.and.ihurrmulti.eq.1) &
            tkerho(i,j,k)=xvuse(i,j,k,nv-1)/volgas(i,j,k)
            if(isemi.eq.1.and.tke(i,j,k).lt.0) tke(i,j,k)=xes(i,j,k,nv-1)
            if(isemi.eq.1.and.tkerho(i,j,k).lt.0) tkerho(i,j,k)=xe(i,j,k,nv-1)
          enddo
        enddo
      enddo

      tauij=0.
      tautilda=0.
      cdrag=0.
      dtemp=0.
      dqw=0.

! compute tautilda(i0,*)

  do i0=1,3    ! i0 denotes x, y, or z direction
  iskip=0
  if(j3.eq.0.and.i0.eq.2) iskip=1
  if(iskip.eq.0) then

! compute tautilda(i0,1) 
  do i=1,np+1
          im1=i-1
          iuse=i
          if(leftedge.eq.1.and.i.eq.1) im1=ipl
          if(rightedge.eq.1.and.i.eq.np+1) iuse=ipr
    do j=1,mp
      ibc=1
      jm1=j-1
      jp1=j+1
          if(topedge.eq.1.and.j.eq.mp) jp1=ipt
          if((topedge.eq.1.and.j.eq.mp).and.isphere.eq.1) ibc=-1
          if(botedge.eq.1.and.j.eq.1)  jm1=ipb
          if((botedge.eq.1.and.j.eq.1).and.isphere.eq.1)  ibc=-1
          if(j3.eq.0) jm1=1
          if(j3.eq.0) jp1=1
      do k=1,l
!     dampx=0.5*(1.+tanh(0.005*(abs(xcrdatas(i,j,k))-(500.*1.e3) )))
!     dampy=0.5*(1.+tanh(0.005*(abs(ycrdatas(i,j,k))-(500.*1.e3) )))
      dampz=0.
!     if(i0.eq.3) &
      dampz=0.5*(1.+tanh(0.005*(zcrdatas(1,1,k)-18000.)))
            ia=(npos-1)*np+i; ja=(mpos-1)*mp+j
!           if(i0.eq.1.and.(ja.eq.170.and.k.eq.65))  &
          kp1=min(k+1,l); km1=max(1,k-1)
          du0(1,:)=dxi*(vvec(iuse,j,k,:)-vvec(im1,j,k,:))
          velyp=0.5*(vvec(iuse,jp1,k,:)+vvec(im1,jp1,k,:))
          velyp(2)=velyp(2)*ibc
          velym=0.5*(vvec(iuse,jm1,k,:)+vvec(im1,jm1,k,:))
          velym(2)=velym(2)*ibc
          du0(2,:)=0.5*dyi*(velyp-velym)

          velzp=0.5*(vvec(iuse,j,kp1,:)+vvec(im1,j,kp1,:))
          velzm=0.5*(vvec(iuse,j,km1,:)+vvec(im1,j,km1,:))
          gii=1./(0.5*(gi(iuse,j,k)+gi(im1,j,k)))

!         if(k.gt.1) then

          du0(3,:)=dzi*(velzp-velzm)/(kp1-km1)

!         else
!         du0(3,:)=dzi*0.5*(velzp+velzm)
!         endif
          xjac0=0.5*(xjac(iuse,j,k,:,:)+xjac(im1,j,k,:,:))
!         rkmave=0.5*(rkm(iuse,j,k)+rkm(im1,j,k))
          if(km_tensor==0)then
            if(itke.eq.1) then

            rkmave=0.09*0.5*(sax(im1,j,k)+sax(iuse,j,k))* &
            sqrt(abs(0.5*(tke(iuse,j,k)+tke(im1,j,k))))+constantdiff

            dxnewa=xjaci(iuse,j,k,1,1)*dx
            dxnewb=xjaci(im1,j,k,1,1)*dx

            rkmave=0.001*0.5*(dxnewa*dxnewa+dxnewb*dxnewb)

!           if(ia.ge.30.and.ia.le.170) then
!           if(ja.ge.30.and.ja.le.170) then
!           rkmave=0.0001*0.5*(dxnewa*dxnewa+dxnewb*dxnewb)
!           endif
!           endif
!           if(i0.eq.3.and.(ia.eq.10.and.ja.eq.10)) print*,k,rkmave, &
!           & 0.0001*0.5*(dxnewa*dxnewa+dxnewb*dxnewb)
            tkerhotmp=0.5*(tkerho(im1,j,k)+tkerho(iuse,j,k))
            if(ismag.eq.1) rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i+1,j,k))+rkmave
!            if(ismag.eq.1) rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i+1,j,k))

            rkmave = 200.0

            rlength=rkmave*0.5*(rhot(iuse,j,k)+rhot(im1,j,k))
            tauij(i,j,k,1)=rlength*tauijf(xjac0,du0,i0,1)
            tautilda(i,j,k,1)=tauijtildatkenew(xjac0,du0,tkerhotmp,rlength,i0,1)
            endif
            if(itke.eq.0) then
            rkmave=0.25*(rkm(iuse,j,k)+rkm(im1,j,k)) &
                      *(rhot(iuse,j,k)+rhot(im1,j,k))
            tauij(i,j,k,1)=rkmave*tauijf(xjac0,du0,i0,1)
            tautilda(i,j,k,1)=rkmave*tauijtilda(xjac0,du0,i0,1)
            endif
          else
! need to load km tensor
!           call getkmtensor(xjac0,dr0,du0,rkmave,kmtensor)
            xsum=0.
            do mm=1,3
              xsum=xsum+kmtensor(i0,mm)*tauijtildatkenew(xjac0,du0,tkerhotmp,rlength,mm,1)
            enddo
            tautilda(i,j,k,1)=xsum
          endif
!       endif
        tautilda(i,j,k,1)=tautilda(i,j,k,1)*gii
      enddo
    enddo
  enddo
!     call mpi_finalize(ierr)
!     stop
! goto 10

! compute tautilda(i0,3)
  do i=1,np
    do j=1,mp
          ibc=1
          if((topedge.eq.1.and.j.eq.mp).and.isphere.eq.1) ibc=-1
          if((botedge.eq.1.and.j.eq.1).and.isphere.eq.1)  ibc=-1
      do k=1,l
        ip1=i+1; im1=i-1
        if(rightedge==1.and.i==np)ip1=ipr
        if(leftedge==1.and.i==1)im1=ipl
        jp1=j+1; jm1=j-1
        if(topedge==1.and.j==mp)jp1=ipt
        if(botedge==1.and.j==1)jm1=ipb
          if(j3.eq.0) jm1=1
          if(j3.eq.0) jp1=1
        km1=max(1,k-1)
          rnoslip=1.
          if(k.eq.1) rnoslip=1.
          velxp=0.5*(vvec(ip1,j,k,:)+vvec(ip1,j,km1,:))*rnoslip
          velxm=0.5*(vvec(im1,j,k,:)+vvec(im1,j,km1,:))*rnoslip
          du0(1,:)=0.5*dxi*(velxp-velxm)
          velyp=0.5*(vvec(i,jp1,k,:)+vvec(i,jp1,km1,:))*rnoslip
          velyp(2)=ibc*velyp(2)
          velym=0.5*(vvec(i,jm1,k,:)+vvec(i,jm1,km1,:))*rnoslip
          velym(2)=ibc*velym(2)
          du0(2,:)=0.5*dyi*(velyp-velym)

          if(k.gt.1) then
          du0(3,:)=dzi*(vvec(i,j,k,:)-vvec(i,j,km1,:))
	  else
	  du0(3,:) = 0.
	  endif

! assumes velocity at ground equal zero...use estimate at middle for ground
!          du0(3,:)=dzi*0.5*(vvec(i,j,2,:)+vvec(i,j,1,:))
!          du0(3,:)=0.		!free-slip

! another approximation...
!         du0(3,:)=dzi*0.5*(vvec(i,j,1,:))
!          du0(3,:)=dzi*(vvec(i,j,2,:)-vvec(i,j,1,:))
!         endif

!	  du0(3,:) = 0.
!	  du0 = 0.

          xjac0=0.5*(xjac(i,j,k,:,:)+xjac(i,j,km1,:,:))
!         rkmave=0.5*(rkm(i,j,k)+rkm(i,j,km1))
!         if(k.eq.1) rkmave=100.0
          if(km_tensor==0)then
            if(itke.eq.1) then

            rkmave=0.09*0.5*(sax(i,j,km1)+sax(i,j,k)) &
                  *sqrt(abs(0.5*(tke(i,j,k)+tke(i,j,km1))))+constantdiff

            dznewa=xjaci(i,j,km1,3,3)*dz
            dznewb=xjaci(i,j,k,3,3)*dz
            rkmave=0.01*0.5*(dznewa*dznewa+dznewb*dznewb)

!           rkmave=0.01*0.5*(rkmax(i,j,k)+rkmax(i,j,km1))
!           rkmave=0.5*(rkmax(i,j,k)+rkmax(i,j,km1))
!           tkerhotmp=0.5*(tkerho(i,j,k)+tkerho(i,j,km1))
!           if(k.eq.1) rkmave=0.075*rkmave+0.1
      if(k.eq.1) windspeed=sqrt((xvuse(i,j,k,1)/xvuse(i,j,k,nv))**2 &
                               +(xvuse(i,j,k,2)/xvuse(i,j,k,nv))**2)
            if(k.eq.1) then
            qv=xvuse(i,j,k,5)/xvuse(i,j,k,nv)
            cpmix=cp*(1.+0.94*qv)
            cvmix=cv*(1.+1.07*qv)
            rgmix=rg*(1.+0.61*qv)
            prrcp=1.e5**(rgmix/cpmix)
            prestmp=(xvuse(i,j,k,4)*rgmix/prrcp)**(cpmix/cvmix)
            pruse=0.5*(1.+tanh((980.*100.-prestmp)/0.1))
            endif

            if(k.eq.1) rkmave=1.*tanh(windspeed/80.)
            if(k.eq.1) damp=0.5*(1.+tanh(0.5*(25-windspeed)))
            if(k.eq.1) rkmave=1.0*(25.*tanh(windspeed/25.)*damp+1.)*pruse &
                             +10.*tanh(windspeed/80)*(1-pruse) 

            if(k.eq.1.and.(irita.eq.1.and.zs(i,j).gt.0.)) & 
            rkmave=250.*tanh(windspeed/80.)

            if(ismag.eq.1.and.k.gt.1) &
            rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i,j,k+1))
!            rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i,j,k+1))+rkmave
!           if(ismag.eq.1.and.k.eq.1) &
!           rkmave=rkmsmag(i,j,k)+rkmave

	    rkmave = 200.0

            rlength=rkmave*0.5*(rhot(i,j,k)+rhot(i,j,km1))
            tauij(i,j,k,3)=rlength*tauijf(xjac0,du0,i0,3)
            tautilda(i,j,k,3)=tauijtildatkenew(xjac0,du0,tkerhotmp,rlength,i0,3)
            endif


            if(itke.eq.0) then
            rkmave=0.25*(rkm(i,j,k)+rkm(i,j,km1))*(rhot(i,j,k)+rhot(i,j,km1))
            if(k.eq.1) rkmave=00.1
!           if(mpi_rank.eq.0.and.i.eq.10) print*,k,rkmave
            tauij(i,j,k,3)=rkmave*tauijf(xjac0,du0,i0,3)
            tautilda(i,j,k,3)=rkmave*tauijtilda(xjac0,du0,i0,3)
!           if((k.eq.38).and.mpi_rank.eq.0)  &
!           print*,i,tautilda(i,j,k,3),rkmave,xjac0,du0
            endif
          else
!           call getkmtensor(xjac0,dr0,du0,rkmave,kmtensor)
            zsum=0.
            do mm=1,3
              zsum=zsum+kmtensor(i0,mm)*tauijtildatkenew(xjac0,du0,tkerhotmp,rlength,mm,3)
            enddo
            tautilda(i,j,k,3)=zsum
          endif
          gii=1./(0.5*(gi(i,j,k)+gi(i,j,km1)))
          tautilda(i,j,k,3)=tautilda(i,j,k,3)*gii
      enddo
    enddo
  enddo

!   	For free slip along lower boundary...
	if (i0 .eq. 3) then
    	do i=1,np
    	do j=1,mp
    	tautilda(i,j,2,3) = tautilda(i,j,3,3) - (tautilda(i,j,4,3) - tautilda(i,j,3,3))
    	tautilda(i,j,1,3) = tautilda(i,j,3,3) - 2.*(tautilda(i,j,4,3) - tautilda(i,j,3,3))
    	enddo
	enddo
	endif



    if(j3.eq.1) then
    do i=1,np
        ip1=i+1; im1=i-1
        if(rightedge==1.and.i==np)ip1=ipr
        if(leftedge==1.and.i==1)im1=ipl
      do j=1,mp+1
          ibc=1
          if((topedge.eq.1.and.j.eq.mp+1).and.isphere.eq.1) ibc=-1
          if((botedge.eq.1.and.j.eq.1).and.isphere.eq.1)  ibc=-1
      juse=j
      jm1=j-1
        if(topedge==1.and.j==mp+1) juse=ipt
        if(botedge==1.and.j==1)jm1=ipb
        do k=1,l
!     dampx=0.5*(1.+tanh(0.005*(abs(xcrdatas(i,j,k))-(500.*1.e3) )))
!     dampy=0.5*(1.+tanh(0.005*(abs(ycrdatas(i,j,k))-(500.*1.e3) )))
      dampz=0.
!     if(i0.eq.3) &
      dampz=0.5*(1.+tanh(0.005*(zcrdatas(1,1,k)-18000.)))
          kp1=min(k+1,l); km1=max(1,k-1)
            velxp(1)=0.5*(vvec(ip1,juse,k,1)+vvec(ip1,jm1,k,1))
            velxp(2)=0.5*(ibc*vvec(ip1,juse,k,2)+ibc*vvec(ip1,jm1,k,2))
            velxp(3)=0.5*(vvec(ip1,juse,k,3)+vvec(ip1,jm1,k,3))
            velxm(1)=0.5*(vvec(im1,juse,k,1)+vvec(im1,jm1,k,1))
            velxm(2)=0.5*(vvec(im1,juse,k,2)+vvec(im1,jm1,k,2))
            velxm(3)=0.5*(vvec(im1,juse,k,3)+vvec(im1,jm1,k,3))

            du0(1,:)=0.5*dxi*(velxp-velxm)
            du0(2,1)=dyi*(vvec(i,juse,k,1)-vvec(i,jm1,k,1))
            du0(2,2)=dyi*(ibc*vvec(i,juse,k,2)-ibc*vvec(i,jm1,k,2))
            du0(2,3)=dyi*(vvec(i,juse,k,3)-vvec(i,jm1,k,3))

            velzp(1)=0.5*(vvec(i,juse,kp1,1)+vvec(i,jm1,kp1,1))
            velzp(2)=0.5*(ibc*vvec(i,juse,kp1,2)+ibc*vvec(i,jm1,kp1,2))
            velzp(3)=0.5*(vvec(i,juse,kp1,3)+vvec(i,jm1,kp1,3))

            velzm(1)=0.5*(vvec(i,juse,km1,1)+vvec(i,jm1,km1,1))
            velzm(2)=0.5*(ibc*vvec(i,juse,km1,2)+ibc*vvec(i,jm1,km1,2))
            velzm(3)=0.5*(vvec(i,juse,km1,3)+vvec(i,jm1,km1,3))

!           if(k.gt.1) then
!!!!            du0(3,:)=0.5*dzi*(velzp-velzm)
            du0(3,:)=dzi*(velzp-velzm)/(kp1-km1)	!new line, above was bug.
!           else
!           du0(3,:)=dzi*0.5*(velzp+velzm)
!           endif
            xjac0=0.5*(xjac(i,juse,k,:,:)+xjac(i,jm1,k,:,:))
!           rkmave=0.5*(rkm(i,juse,k)+rkm(i,jm1,k))
            if(km_tensor==0)then
            if(itke.eq.1) then

            rkmave=0.09*0.5*(sax(i,jm1,k)+sax(i,juse,k)) &
                  *sqrt(abs(0.5*(tke(i,juse,k)+tke(i,jm1,k))))+constantdiff
            dynewa=xjaci(i,juse,k,2,2)*dy 
            dynewb=xjaci(i,jm1,k,2,2)*dy 
            ia=(npos-1)*np+i; ja=(mpos-1)*mp+j

            rkmave=0.001*0.5*(dynewa*dynewa+dynewb*dynewb)

!           if(ia.ge.30.and.ia.le.170) then
!           if(ja.ge.30.and.ja.le.170) then
!           rkmave=0.0001*0.5*(dynewa*dynewa+dynewb*dynewb)
!           endif
!           endif
!           if(i0.eq.3.and.(ia.eq.10.and.ja.eq.10)) print*,k,rkmave, &
!           & 0.0001*0.5*(dxnewa*dxnewa+dxnewb*dxnewb)
            tkerhotmp=0.5*(tkerho(i,jm1,k)+tkerho(i,juse,k))
            if(ismag.eq.1)  rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i,j+1,k))+rkmave
!            if(ismag.eq.1)  rkmave=0.5*(rkmsmag(i,j,k)+rkmsmag(i,j+1,k))

	    rkmave = 200.0

            rlength=rkmave*0.5*(rhot(i,juse,k)+rhot(i,jm1,k))
            tauij(i,j,k,2)=rlength*tauijf(xjac0,du0,i0,2)
            tautilda(i,j,k,2)=tauijtildatkenew(xjac0,du0,tkerhotmp,rlength,i0,2)
            endif
            if(itke.eq.0) then
            rkmave=0.25*(rkm(i,juse,k)+rkm(i,jm1,k))*(rhot(i,juse,k)+rhot(i,jm1,k))
            tauij(i,j,k,2)=rkmave*tauijf(xjac0,du0,i0,2)
              tautilda(i,j,k,2)=rkmave*tauijtilda(xjac0,du0,i0,2)
            endif
            else
!             call getkmtensor(xjac0,dr0,du0,rkmave,kmtensor)
              ysum=0.
              do mm=1,3
                ysum=ysum+kmtensor(i0,mm)*tauijtildatkenew(xjac0,du0,tkerho,rlength,mm,2)
              enddo
              tautilda(i,j,k,2)=ysum
            endif
!         endif
          gii=1./(0.5*(gi(i,juse,k)+gi(i,jm1,k)))
          tautilda(i,j,k,2)=tautilda(i,j,k,2)*gii
        enddo
      enddo
    enddo
    endif

! compute shear generation terms
     if(isweep.eq.1.and.itke.eq.1) &
     call tkegennew(vvec(1-ih,1-ih,1,i0),tauij,1-ih,np+ih,1-ih,mp+ih,l,i0)

    do i=1,np
      do j=1,mp
        do k=1,l
          kp0=min(k+1,l); km0=k
          ip0=i+1; im0=i
          jp0=j+1; jm0=j
          if(j3.eq.0) jp0=1
          txp=tautilda(ip0,j,k,1)
          txm=tautilda(im0,j,k,1)
          typ=tautilda(i,jp0,k,2)
          tym=tautilda(i,jm0,k,2)
          tzp=tautilda(i,j,kp0,3)
          tzm=tautilda(i,j,km0,3)
!          if(k==l)tzp=0.
          tx=dxi*(txp-txm)
          ty=dyi*(typ-tym)
          tz=dzi*(tzp-tzm)
          ftau=(tx+ty+tz)*gi(i,j,k)

            if(isweep.eq.1) &
            fdissip(i,j,k,i0)=-ftau*dt

            ia=(npos-1)*np+i; ja=(mpos-1)*mp+j

            if(isweep.eq.2) &
            fdissip(i,j,k,i0+4)=-ftau*dt


!	if (time .eq. 20 .and. i0.eq.3.and.ia.eq.170.and.ja.eq.170.and.k.le.4)then
!	print*,k,fdissip(i,j,k,i0),tx,ty,tz
!	endif

        enddo
      enddo
    enddo

!	if(time .eq. 20 .and. i0.eq.3) then
!	call mpi_finalize(ierr)
!	stop
!	endif

!    	do i=1,np
!  	do j=1,mp
!	fdissip(i,j,1:20,3) = fdissip(i,j,20,3)
!	enddo
!	enddo



! endif
! enddo

 endif
 enddo   ! i0 loop
! call rmaxmin1(fdissip(1-ih,1-ih,1,1),'f1',1-ih,np+ih,1-ih,mp+ih,l)
! call rmaxmin1(fdissip(1-ih,1-ih,1,2),'f2',1-ih,np+ih,1-ih,mp+ih,l)
! call rmaxmin1(fdissip(1-ih,1-ih,1,3),'f3',1-ih,np+ih,1-ih,mp+ih,l)
!     call mpi_finalize(ierr)
!     stop
100   continue
!     call mpi_finalize(ierr)
!     stop

! if(made_tmp_zero)deallocate(zero)


!  if(time .eq. 20) then
!  do kv=1,3
!  call writeio(fdissip(1-ih,1-ih,1,kv),'u',1,new,21,fname,1-ih,np+ih,1-ih,mp+ih,l)
!  enddo

!  do kv=1,3
!  call writeio(tautilda(1-ih,1-ih,1,kv),'u',1,new,21,fname,1-ih,np+ih,1-ih,mp+ih,l)
!  enddo

!	call mpi_finalize(ierr)
!	stop
!  endif

  deallocate(tautilda,tauij,tke,tkerho)
  deallocate (vvec)
  deallocate (rhot)

!	if(mpi_rank.eq.0)then
!	print*,fdissip(100,100,:,1)
!	endif

!	call mpi_finalize(ierr)
!	stop

!        if (mpi_rank.eq.0)then
!        print*,'Dissipation values in height:'
!        print*,fdissip(10,10,:,1)
!        endif


  end