--- From the main body of triple_zyx.f -------------------------------
***
      call xopause(tbar,beta,zip,work,work(nz+2),nz,ibz,rx,ry,rz,ipro)
***
          call wavepacket(tt,ww,v3,tbar,beta,work,work(nxp+1),
     &                    work(nxp+nz+2),nx,ny,nz,n3,ibx,iby,ibz,
     &                 tgx,tgy,tgz,wp,fn,mp,ipro,mype,ncpu,ncput)
***
 --- And the subroutine -----------------------------------------------
!***********************************************************************
      subroutine xopause(tbar,beta,zip,z,a,nz,ibz,rx,ry,rz,iprob)
!-----------------------------------------------------------------------
!  Initialize background T in physical space for an x-o-pause.
!
!  If more than one of rx,ry,rz is non-zero, i.e., the box is tilted, 
!  then bail with an error message.
!
!  Actually, we only allow rz=1, rx=ry=0.
!
!  The background temperature we set here is:
!
!  dT/dz = 1.0 + (tzr-1)*.5*{
!                  tanh[(z-zmx1)*.5/t]  -  tanh[(z-zmx)*.5/t]
!                + tanh[(z+zmx2)*.5/t]  -  tanh[(z    )*.5/t]  }
!
!  This corresponds to the temperature profile:
!
!  T = z + (tzr-1)*t*{
!            ln{cosh[(z-zmx1)*.5/t]}  -  ln{cosh[(z-zmx)*.5/t]}
!          + ln{cosh[(z+zmx2)*.5/t]}  -  ln{cosh[(z    )*.5/t]}  }
!
!  Note: integral{tanh((z-z0)/2t)} = 2t*ln{cosh[(z-z0)/2t]}
!
! _____ _____ _____________________________________________ top
!   |     |    |                               .^`
!   |     |    |                            .^`
!   |     |    |                         .^`         
!   |    zmx2  |                      .^`   <== dT/dz = tzr
!   |     |    |                   .^`         
!   |     |    |                .^`        
!   |   __|__  |_____________.^`___________________________ x-o-pause
!   |     |    |            /
!   |     |    |           /
!  zmx    |    |          /
!   |     |    |         /
!   |     |    |        /
!   |    zmx1  |       /
!   |     |    |      /  <== normalized temperaure gradient, dT/dz = 1
!   |     |    |     /
!   |     |    |    /
!   |     |    |   /
!   |     |    |  /
!   |     |    | /
! __|__ __|__ _|/__________________________________________ bottom
!
!
!             beta  =  [T(top)-T(bottom)] / zmx
!
!-----------------------------------------------------------------------
      ! include some parameter files
      include 'include/garten.fi'
      include 'include/pars.fi'

      ! declare arrays
      dimension z(nz+1),a(nz+1,4),tbar(nz+ibz,2),zip(nz+ibz)

      nzt = nz + ibz
      do k=1,nzt
        zip(k) = 0.0
      enddo

      if ((iprob.ne.12).and.(iprob.ne.13)) then
        do k=1,nzt
          tbar(k,1) = 0.0
          tbar(k,2) = 1.0
        enddo
        beta = 1.0
        return
      endif

      ! check that only rz is non-zero
      if ((rx.ne.0).or.(ry.ne.0).or.(rz.eq.0)) then
        write(6,*)"  * * * * * * * * * * * * * * * * * * * * * * * * *"
        write(6,*)"  *                                               *"
        write(6,*)"  *       you called subroutine xopause with      *"
        write(6,*)"  *         rx.ne.0 or ry.ne.0 or rz.eq.0         *"
        write(6,*)"  *                                               *"
        write(6,*)"  *                BAILING OUT!                   *"
        write(6,*)"  *                                               *"
        write(6,*)"  * * * * * * * * * * * * * * * * * * * * * * * * *"
        stop
      endif

      ! define some useful constants
      zmx1 = zmx - zmx2
      f1   = (tzr-1.)*.5

      ! define z(nzt)
      do k=1,nz+1
        z(k) = float(k-1)/float(nz)*zmx
      enddo
      delz = z(nz+1) - z(1)
     
      ! obtain tanh profiles
      call tanhprof(a(1,1),z, zmx1,edge,nz+1)
      call tanhprof(a(1,2),z, zmx ,edge,nz+1)
      call tanhprof(a(1,3),z,-zmx2,edge,nz+1)
      call tanhprof(a(1,4),z, 0.0 ,edge,nz+1)
      
      ! build mean dtdz profile
      do k=1,nzt
        tbar(k,2) = 1.0  + f1*( a(k,1) - a(k,2) + a(k,3) - a(k,4) )
      enddo

      ! obtain integrals of tanh profiles
      call tanhint(a(1,1),z, zmx1,edge,nz+1)
      call tanhint(a(1,2),z, zmx ,edge,nz+1)
      call tanhint(a(1,3),z,-zmx2,edge,nz+1)
      call tanhint(a(1,4),z, 0.0 ,edge,nz+1)
      
      ! define modified limit values for a(k,1)-a(k,2)+a(k,3)-a(k,4)
      ap    = -f1*( a(nz+1,1) - a(nz+1,2) + a(nz+1,3) - a(nz+1,4) )/delz
      am    = -f1*( a(   1,1) - a(   1,2) + a(   1,3) - a(   1,4) )/delz
      delaz = am*z(nz+1) - ap*z(1)
      dela  = ap - am

      ! the mean temperature gradient from bottom to top;
      ! use this instead of 1.
      beta  = 1. - dela

      ! build mean T profile
      do k=1,nzt
        tbar(k,1) = delaz + dela*z(k) + f1*(a(k,1)-a(k,2)+a(k,3)-a(k,4))
      enddo
!-----------------------------------------------------------------------
      return
      end


 --- Here are the 2 other subroutines used in computing tbar -----------
!***********************************************************************
      subroutine tanhint(a,z,z1,t,nzt)
!-----------------------------------------------------------------------
!  Compute the analytic integral of tanh((z-z1)/t)
!-----------------------------------------------------------------------
      dimension a(nzt),z(nzt)
!-----------------------------------------------------------------------
      ti = 1./t
      t2 = 2.*t
      do k=1,nzt
        zdiff = z(k)-z1
        if (zdiff.lt.0.) then
          a(k)=-zdiff+t2*log(1+exp( zdiff*ti))
        else
          a(k)= zdiff+t2*log(1+exp(-zdiff*ti))
        endif
      enddo
!-----------------------------------------------------------------------
      return
      end
!***********************************************************************
      subroutine tanhprof(da,z,z1,t,nzt)
!-----------------------------------------------------------------------
!  Return tanh((z-z1)/t)
!-----------------------------------------------------------------------
      dimension da(nzt),z(nzt)
!-----------------------------------------------------------------------
      t2i = .5/t
      do k=1,nzt
        zdiff = z(k)-z1
        da(k)=tanh(zdiff*t2i)
      enddo
      return
      end