!******************** start file outptb.for ; group outptb *************
!  outptb -- renormalize beam mc code outputs after orbit loop
!
!  dmc-- reorganized 30 aug 1988, in support of conversion of code to
!  multispecies beams.
!
!  some code formerly located in this routine has been moved to outptb2
!
!    dmc july 1998 -- to support higher resolution radial grid in rest
!    of code, smoothing of radial profiles is now supported in the mc
!    code.  ** profiles renormalized in outptb should also be smoothed
!    in subroutine orbsmoo **
!    see also:  outptb2, orbsmoo2.
!    and see also:  depnrm, depsmoo.
!
subroutine outptb(iorb,ierr)
  !
  !  iorb = no. of orbits traced; if iorb=0, no orbits were followed.
  !
  use nbi_com
  use nbi_mpi_com
  !
  !  *** modif dmc 5 mar 1985 ***
  !   as preparatory step for generalized geometry beam code,
  !   most outputs are having their volume normalization removed, because
  !   geometry (hence volume) can change from beam timestep to beam
  !   timestep.  therefore outputs have units such as n, ev, watts, n/sec,
  !   **not** n/cm3, j/cm3, watts/cm3, n/cm3/sec
  !
  !  --> as of 5 mar, the poloidal profile outputs and the beam dist. fcn
  !      still have their volume normalization
  !
  !
  !	this subroutine outputs source terms
  !	from the monte-carlo beam orbit code
  !
  !		upgrade to distinguish cx loss from inside the
  !		plasma from cx loss from outside (bpcxin, bpcxot
  !		rather than bpcx)  rob goldston july 21, 1983
  !
  !  constants used in normalization
  !
  !
  !============
  implicit none
  integer :: iorb,j,jp1,ig,ien,ie,i,ip1,k,ierr
  !============
  real(fp) :: zdti,zbpf,zdti7,zomega,ztircor,znvz,znv,zjldot
  real(fp) :: zbmnrm,zbmn
  real(fp) :: zdelebd,zvp,zdp,zvpp,zdpp
  !============

  real(fp) :: fbm_flr   ! fbm selector function

  real(fp) :: bsumint,zeav,zeav2,zde2,znt
  integer :: mk,ib,ii
  !
  !-----------------------------------------------------------------
  !
  ierr=0
  !
  !  timestep normalization
  !
  zdti=1.0_fp/delt
  !
  !  turn ptcl*ev into joules
  !
  zbpf=xe*zdti
  !
  !  normalize ang. momementum balance terms
  !
  zdti7=1.0e-7_fp*zdti

  !
  !  normalize beam scalar quantities
  !
  !-----------------------------
  !  following may also have contributions from dep/depnrm
  !
  !
  bpls(lsbeam)=(bpls(lsbeam)+bpls_dep(lsbeam))*zbpf ! orbit loss power...
  bpls_dep(lsbeam)=bpls_dep(lsbeam)*zbpf
  !
  bxflws(lsbeam)=(bxflws(lsbeam)+bxflws_dep(lsbeam))*zdti ! ptcls...
  bxflws_dep(lsbeam)=bxflws_dep(lsbeam)*zdti
  bxflws_prompt(lsbeam)=bxflws_prompt(lsbeam)*zdti+ &
       bxflws_dep(lsbeam)

  bphorbs(lsbeam)=(bphorbs(lsbeam)+bphorbs_dep(lsbeam))*zdti7 ! momentum
  bphorbs_dep(lsbeam)=bphorbs_dep(lsbeam)*zdti7
  !
  bphdeps(lsbeam)=bphdeps(lsbeam)*zdti7 ! mo. deposition at ptcl position

  bphdepsgc(lsbeam)=bphdepsgc(lsbeam)*zdti7 ! mo. deposition at gc

  bphinjs(lsbeam)=bphinjs(lsbeam)*zdti7 ! mo. injected

  bphshins(lsbeam)=bphshins(lsbeam)*zdti7 ! shine through mo.

  !
  bphjxbs(lsbeam)=(bphjxbs(lsbeam)+bphjxbs_ntr(lsbeam)+ &
       bphjxbs_dep(lsbeam))*zdti7 ! jxb torque
  bphjxbs_dep(lsbeam)=bphjxbs_dep(lsbeam)*zdti7
  !
  bpjxbs(lsbeam)=(bpjxbs(lsbeam)+bpjxbs_ntr(lsbeam)+ &
       bpjxbs_dep(lsbeam))*zbpf ! jxb work
  bpjxbs_dep(lsbeam)=bpjxbs_dep(lsbeam)*zbpf

  ! momentum balance check
  bphdpba(lsbeam)=bphinjs(lsbeam)-bphdeps(lsbeam)-bphshins(lsbeam)

  pbdepmc(lsbeam)=pbdepmc(lsbeam)*zbpf

  ! power deposited balance check
  pbdepba(lsbeam)=pinjs(lsbeam)-pbdepmc(lsbeam)-bpshins(lsbeam)


  !
  !-----------------------------
  !
  !  following scalars are summed in orbit loop
  !
  !  fishbone losses
  bpfshs(lsbeam)=bpfshs(lsbeam)*zbpf
  bxfshs(lsbeam)=bxfshs(lsbeam)*zdti
  bphfshs(lsbeam)=bphfshs(lsbeam)*zdti7
  !  fishbone lossesaep (mhd activity) power

  pbaepms(1:mj,1:mnkaep,lsbeam)=pbaepms(1:mj,1:mnkaep,lsbeam)*zbpf
  !
  !  ripple losses
  bprpls(lsbeam)=bprpls(lsbeam)*zbpf
  bxrpls(lsbeam)=bxrpls(lsbeam)*zdti
  bphrpls(lsbeam)=bphrpls(lsbeam)*zdti7
  !
  !  charge exchange losses
  snbcxs(lsbeam)=(snbxtot(lsbeam)-snbcxos(lsbeam))*zdti
  snbcxos(lsbeam)=snbcxos(lsbeam)*zdti  ! from outside plasma
  sbcxxs(lsbeam)=sbcxxs(lsbeam)*zdti    ! from inside plasma
  !
  bpcxins(lsbeam)=bpcxins(lsbeam)*zbpf
  bpcxots(lsbeam)=bpcxots(lsbeam)*zbpf
  bpcxi0s(lsbeam)=bpcxi0s(lsbeam)*zbpf
  bpcxx0s(lsbeam)=bpcxx0s(lsbeam)*zbpf
  bpcxers(lsbeam)=bpcxers(lsbeam)*zbpf
  bpcxris(lsbeam)=bpcxris(lsbeam)*zbpf
  bpcxrxs(lsbeam)=bpcxrxs(lsbeam)*zbpf
  !
  !  cx momentum losses
  bph0was(lsbeam)=bph0was(lsbeam)*zdti7
  bphcxls(lsbeam)=bphcxls(lsbeam)*zdti7
  bphrcas(lsbeam)=bphrcas(lsbeam)*zdti7
  !
  !  momentum terms
  bphcols(lsbeam)=(bsumint(tqbes(1:mj,lsbeam))+ &
       bsumint(tqbis(1:mj,lsbeam)))*zdti7
  bphths(lsbeam)=bsumint(tqbths(1:mj,lsbeam))*zdti7
  bphohbs(lsbeam)=bsumint(tqohbs(1:mj,lsbeam))*zdti7
  !
  !  orbit code error checks
  bpherrs(lsbeam)=bpherrs(lsbeam)*zdti7
  bperrs(lsbeam)=bperrs(lsbeam)*zbpf

  !
  bpephis(lsbeam)=0.0_fp
  !
  !  anomalous diffusion momentum terms
  bphdfbs(lsbeam)=bphdfbs(lsbeam)*zdti7
  !
  bptdfbs(lsbeam)=bptdfbs(lsbeam)*zbpf
  !
  !  rotating source friction term
  bpthsfs(lsbeam)=bpthsfs(lsbeam)*zbpf

  !
  !----------------------------------------
  !  normalize profile quantities, and add up some scalar quantities
  !
  bpths(lsbeam)=0.0_fp
  bpohs(lsbeam)=0.0_fp
  !
  bptcols(lsbeam)=0.0_fp
  bpthros(lsbeam)=0.0_fp
  bpthafs(lsbeam)=0.0_fp
  !
  do j=lcentr,ledge
     jp1=j+1
     !
     !  dv*<j dot b> -> da*<j dot b>
     !    will be converted to da*jt form (amps) in outptb2...
     !
     !  all driven currents
     pljbbms(j,lsbeam)=pljbbms(j,lsbeam)*xe*zdti* &
          dareazc_r8(j)/dvolzc_r8(j)
     !
     !  driven from "co" beams only; no fusion product contributions
     pljbcct(j,lsbeam)=pljbcct(j,lsbeam)*xe*zdti* &
          dareazc_r8(j)/dvolzc_r8(j)
     !
     pbis(j,lsbeam)=pbis(j,lsbeam)*zbpf
     pbes(j,lsbeam)=pbes(j,lsbeam)*zbpf
     !
     pbcprss(j,lsbeam)=pbcprss(j,lsbeam)*zbpf
     !
     pbdfbis(j,lsbeam)=pbdfbis(j,lsbeam)*zbpf
     pbdfbes(j,lsbeam)=pbdfbes(j,lsbeam)*zbpf
     !
     pbths(j,lsbeam)=pbths(j,lsbeam)*zbpf
     bpths(lsbeam)=bpths(lsbeam)+pbths(j,lsbeam)
     !
     if(nlbfpp) then
        pbthr2s(j,lsbeam)=0.0_fp  ! dmc -- bugfix 8 jun 2000
     else
        zomega=0.5_fp*(omegag(j)+omegag(j+1))
        pbthr2s(j,lsbeam)=v2toev*zomega**2* &
             (pbthr2s(j,lsbeam)-gr2zc(j)*sorbths(j,lsbeam))*zbpf
     end if
     !
     bpthafs(lsbeam)=bpthafs(lsbeam)+pbthr2s(j,lsbeam)
     !
     pbohs(j,lsbeam)=pbohs(j,lsbeam)*zbpf
     bpohs(lsbeam)=bpohs(lsbeam)+pbohs(j,lsbeam)
     !
     !  don't neglect dep contribution
     sorbes(j,lsbeam)=(sorbes(j,lsbeam)+sorbes_ntr(j,lsbeam)+ &
          sorbes_dep(j,lsbeam))*zdti
     sorbes_dep(j,lsbeam)=sorbes_dep(j,lsbeam)*zdti
     !
     sorbths(j,lsbeam)=sorbths(j,lsbeam)*zdti
     !
     !  nutrav outputs
     sbcx0s(j,lsbeam)=sbcx0s(j,lsbeam)*zdti
     sbcxris(j,lsbeam)=sbcxris(j,lsbeam)*zdti
     sbcxres(j,lsbeam)=sbcxres(j,lsbeam)*zdti
     sbcxrzs(j,lsbeam)=sbcxrzs(j,lsbeam)*zdti
     sbcxrls(j,lsbeam)=sbcxrls(j,lsbeam)*zdti
     sbcxrbs(j,lsbeam)=sbcxrbs(j,lsbeam)*zdti
     !
     do ig=1,ng
        sbcxrxs(ig,j,lsbeam)=sbcxrxs(ig,j,lsbeam)*zdti
        !  dep contribution to be included...
        sorbhs(j,ig,lsbeam)= &
             (sorbhs(j,ig,lsbeam)+sorbhs_ntr(j,ig,lsbeam)+ &
             sorbhs_dep(j,ig,lsbeam))*zdti
        sorbhs_dep(j,ig,lsbeam)=sorbhs_dep(j,ig,lsbeam)*zdti
        sorbhs_ntr(j,ig,lsbeam)=sorbhs_ntr(j,ig,lsbeam)*zdti
        !
        sorbn0s(j,ig,lsbeam)= &
             (sorbn0s(j,ig,lsbeam)+sorbn0s_dep(j,ig,lsbeam))*zdti
        sorbn0s_dep(j,ig,lsbeam)=sorbn0s_dep(j,ig,lsbeam)*zdti
        !
        sinbn0s(j,ig,lsbeam)=sinbn0s(j,ig,lsbeam)*zdti
        sinb0is(j,ig,lsbeam)=sinb0is(j,ig,lsbeam)*zdti


     end do
     !
     !------------------
     !  torques
     !
     tqohbs(j,lsbeam)=tqohbs(j,lsbeam)*zdti7
     !
     if(wtqbxfs(j,lsbeam).ne.0.0_fp) then
        tqbxfs(j,lsbeam)=tqbxfs(j,lsbeam)/wtqbxfs(j,lsbeam)
     else
        tqbxfs(j,lsbeam)=0.0_fp
     end if
     !
     zomega=0.5_fp*(omegag(j)+omegag(j+1))
     tqbths(j,lsbeam)=tqbths(j,lsbeam)*zdti7
     bpthros(lsbeam)=bpthros(lsbeam)+tqbths(j,lsbeam)*zomega
     !
     tqbis(j,lsbeam)=tqbis(j,lsbeam)*zdti7
     tqbes(j,lsbeam)=tqbes(j,lsbeam)*zdti7
     bptcols(lsbeam)=bptcols(lsbeam)+ &
          (tqbes(j,lsbeam)+tqbis(j,lsbeam))*zomega
     !
     tqdifbs(j,lsbeam)=tqdifbs(j,lsbeam)*zdti7
     !
     !  include compreesion, dep contribution, if any ...
     !
     tqjxbs(j,lsbeam)=(tqjxbs(j,lsbeam)+tqjxbs_dep(j,lsbeam)+ &
          tqjxbs_ntr(j,lsbeam) )*zdti7
     tqjxbs_dep(j,lsbeam)=tqjxbs_dep(j,lsbeam)*zdti7
     !
     tqrpls(j,lsbeam)=tqrpls(j,lsbeam)*zdti7
     !
     tqbcxs(j,lsbeam)=(tqbcxs(j,lsbeam)+tqbcxs_dep(j,lsbeam))*zdti7
     tqbcxs_dep(j,lsbeam)=tqbcxs_dep(j,lsbeam)*zdti7
     !
     pbephis(j,lsbeam)=(pbephis(j,lsbeam)+pbephis_dep(j,lsbeam)+ &
          pbephis_ntr(j,lsbeam) )*zbpf
     pbephis_dep(j,lsbeam)=pbephis_dep(j,lsbeam)*zbpf
     bpephis(lsbeam)=bpephis(lsbeam)+pbephis(j,lsbeam)
     !
     !  radial current sums
     tirtots(jp1,lsbeam)= &
          (tirtots(jp1,lsbeam)+tirtots_dep(jp1,lsbeam)+ &
          tirtots_ntr(jp1,lsbeam) )*zbpf
     tirtots_dep(jp1,lsbeam)=tirtots_dep(jp1,lsbeam)*zbpf
     fpixbts(jp1,lsbeam)= &
          (fpixbts(jp1,lsbeam)+fpixbts_dep(jp1,lsbeam)+ &
          fpixbts_ntr(jp1,lsbeam) )*1.e-2_fp*zbpf
     fpixbts_dep(jp1,lsbeam)=fpixbts_dep(jp1,lsbeam)*1.e-2_fp*zbpf
     !
     tirorbs(jp1,lsbeam)=tirorbs(jp1,lsbeam)*zbpf
     tirabds(jp1,lsbeam)=tirabds(jp1,lsbeam)*zbpf
     tirfshs(jp1,lsbeam)=tirfshs(jp1,lsbeam)*zbpf
     tirrips(jp1,lsbeam)=tirrips(jp1,lsbeam)*zbpf
     !
     !  see tirorbs correction, separate loop, below...
     !
     !---------------------------------------------------------------------
     !
     !	define beam avg energy (ions) in ev, and handle beam density
     !	and beta
     !
     if(bdenss(j,lsbeam).gt.0._fp) then
        !  average total, parallel, and perp beam ion energies
        omegabs(j,lsbeam)=omegabs(j,lsbeam)/bdenss(j,lsbeam)
        udenspp(j,lsbeam)=ebapp2s(j,lsbeam)*zbpf
        udenspl(j,lsbeam)=ebapl2s(j,lsbeam)*zbpf
        ebavps(j,lsbeam)=(ebapp2s(j,lsbeam)+ebapl2s(j,lsbeam))/ &
             bdenss(j,lsbeam)
        ebapp2s(j,lsbeam)=ebapp2s(j,lsbeam)/bdenss(j,lsbeam)
        ebapl2s(j,lsbeam)=ebapl2s(j,lsbeam)/bdenss(j,lsbeam)
        xzbmav(j,lsbeam)=xzbmav(j,lsbeam)/bdenss(j,lsbeam)
        bdenss(j,lsbeam)=bdenss(j,lsbeam)*zdti
     else
        omegabs(j,lsbeam)=0.0_fp
        udenspp(j,lsbeam)=0.0_fp
        udenspl(j,lsbeam)=0.0_fp
        ebavps(j,lsbeam)=0.0_fp
        ebapp2s(j,lsbeam)=0.0_fp
        ebapl2s(j,lsbeam)=0.0_fp
        xzbmav(j,lsbeam)=0.0_fp
        bdenss(j,lsbeam)=0.0_fp
     end if
     mc_dens(j,lsbeam) = mc_dens(j,lsbeam)/dvolzc(j)
     !
     !  energy-fraction-resolved outputs
     !
     do ien=1,3
        znt=bnbs_ien(j,ien,lsbeam)
        if(znt.gt.0.0_fp) then
           upbs_ien(j,ien,lsbeam)=upbs_ien(j,ien,lsbeam)/znt
           bnbs_ien(j,ien,lsbeam)=bnbs_ien(j,ien,lsbeam)*zdti
           pbes_ien(j,ien,lsbeam)=pbes_ien(j,ien,lsbeam)*zbpf
           pbis_ien(j,ien,lsbeam)=pbis_ien(j,ien,lsbeam)*zbpf
        end if
     end do
     !
     !  dmc jul 1992
     !
     do ie=1,nerngfi
        densfi(j,ie,lsbeam)=densfi(j,ie,lsbeam)*zdti
     end do
     !
     !  dmc mar 2007 fnb, fvdotb at gc
     !
     ! commented by mg jan 30, 2012
     !        do ie=1,nznbmea(lsbeam)
     !           if(fnb(ie,j,lsbeam).eq.0.0_fp) then
     !              fvdotb(ie,j,lsbeam)=0.0_fp
     !           else
     !              fvdotb(ie,j,lsbeam)=fvdotb(ie,j,lsbeam)/fnb(ie,j,lsbeam)
     !              fnb(ie,j,lsbeam) = fnb(ie,j,lsbeam)*zdti
     !           end if
     !        end do
     !
     !  end of 1d radial profile renorm (j) loop

  end do
  !
  !
  bdenss(lep1,lsbeam)=bdenss(lep1,lsbeam)*zdti
  !
  !  loop to correct the radial current due to orbiting...
  !
  do j=lcentr,ledge
     jp1=j+1
     ztircor=0.0_fp
     !
     !  tirorbs sum records motion of bounce average orbit location;
     !  it may be necessary to subtract off non-orbit effects on the
     !  orbit location e.g. anom diffusion which has its own current
     !  profile already computed...
     !
     !  anom diffusion
     !
     ztircor=ztircor+tirabds(jp1,lsbeam)
     !
     !  xi expansion/contraction relative to flux surfaces
     !  dmc may 2010 -- removed this code, it is not needed now that a
     !  compression operator is run at the start of orbiting; this should
     !  have been removed a long time ago... lines commented with "cdmc"
     !
     !  the correction is d/d(xi) of the volume integral of the charge
     !  density, times xi*eldot, z*n*dv*xi*eldot; the dv has not yet been
     !  divided out of bdenss.  for heavy ions, z is a profile.
     !
     !dmc        if(nlheavy(lsbeam)) then
     !dmc          znvz=0.5_fp*(bdenss(j,lsbeam)*xzbmav(j,lsbeam)+
     !dmc     >                bdenss(jp1,lsbeam)*xzbmav(jp1,lsbeam))
     !dmc        else
     !dmc          znv=0.5_fp*(bdenss(j,lsbeam)+bdenss(jp1,lsbeam))
     !dmc          znvz=znv*xzbeami             ! n*dv*z+
     !dmc        end if
     !dmc        zjldot=xe*xi(jp1,1)*eldot*znvz
     !
     !dmc        ztircor=ztircor+zjldot
     !
     tirorbs(jp1,lsbeam)=tirorbs(jp1,lsbeam)-ztircor
     !
  end do
  !
  !----------------------------------------------------------
  !  abdifu net diffusivity and velocity profiles
  !
  do j=lcentr,ledge
     if(wgvdx(j,lsbeam).gt.zero) then
        ! zonal average db & vb from abdifu
        difbx(j,lsbeam) = difbx(j,lsbeam)/wgvdx(j,lsbeam)
        velbx(j,lsbeam) = velbx(j,lsbeam)/wgvdx(j,lsbeam)
     end if
  end do
  !
  !  shift to zone boundary orientation
  !
  difbx(lep1,lsbeam)=difbx(ledge,lsbeam)
  velbx(lep1,lsbeam)=velbx(ledge,lsbeam)

  zdp=difbx(lcentr,lsbeam)
  zvp=velbx(lcentr,lsbeam)

  do j=lcentr,ledge
     zdpp=zdp
     zvpp=zvp

     zdp=difbx(j,lsbeam)
     zvp=velbx(j,lsbeam)

     difbx(j,lsbeam) = 0.5_fp*(zdp+zdpp)
     velbx(j,lsbeam) = 0.5_fp*(zvp+zvpp)
  end do
  !
  !----------------------------------------------------------
  !
  !  normalize beam distribution		(r.j.g. : 4-11-82)
  !
  nfbmpt(1:mimxbz,lsbeam)=0 ! count non-zero v-space bins per spatial bin
  !
  if(.not.nlfbm .or. nlfbmfpp(lsbeam) .or. nlheavy(lsbeam)) &
       goto 210

  if(iorb.eq.0) goto 210
  !
  !	** oct 1982 ** flip the index organization of fbm for efficient
  !	access for use (as opposed to efficient access for calculation)
  !
  !ccc      call flipr    ! dmc feb 2009 -- trying direct accumulation of
  !                       ! fbm & fbm_ptcl -- therefore, no flipr
  !
  !
  !
  !calculation of bdens2 eba2pl eba2pp vbtor2 -moved to outptb2
  zbmnrm=real(nznbma)/(delt)
  do k=1, nfbznsi
     zbmn=zbmnrm/bmvol(k)
     do j=1, nznbma
        do i=1, nznbmea(lsbeam)
           ip1=i+1
           zdelebd=efbmb(ip1,lsbeam)-efbmb(i,lsbeam)
           fbm(i,j,k,lsbeam) = fbm(i,j,k,lsbeam)*zbmn &
                /zdelebd
           !     at the particle position
           if(nlfbmflr)then
              do ig = 1, mgyro ! gyro phase index
                 fbm_ptcl(i,j,k,lsbeam,ig) = fbm_ptcl(i,j,k, &
                      lsbeam,ig)*zbmn/zdelebd
              end do
           end if !nlfbmflr
           !
           ! fbm_flr
           if(fbm_flr(i,j,k,lsbeam,0).gt.0.0_fp) &
                nfbmpt(k,lsbeam)=nfbmpt(k,lsbeam) + 1
        end do
     end do
  end do
210 continue
  !
  !-----------------------------------------------------------------------
  !  dmc feb 1992 -- addl particle balance diagnostic outputs normalized
  !    here
  !
  !  cx fast neutral sce/recapture
  sfcx0mc(lsbeam)=sfcx0mc(lsbeam)*zdti
  sfcxmc1d(:,lsbeam)=sfcxmc1d(:,lsbeam)*zdti
  sfcxesc(lsbeam)=sfcxesc(lsbeam)*zdti
  sfcxrr(lsbeam)=sfcxrr(lsbeam)*zdti
  sfcxrmc(lsbeam)=sfcxrmc(lsbeam)*zdti
  !
  !  cx fast neutral particle balance
  sfcxbal(lsbeam)=sfcx0mc(lsbeam) &
       -sfcxesc(lsbeam)-sfcxrr(lsbeam)-sfcxrmc(lsbeam)
  !
  !  orbit code particle balance
  sforbal(lsbeam)=sfdepmc(lsbeam)+sfdeprr(lsbeam) &
       -bxflws(lsbeam)-bsumint(sorbths(1:mj,lsbeam)) &
       -bxfshs(lsbeam) &
       -sfcxesc(lsbeam)-sfcxrr(lsbeam)-bxrpls(lsbeam) &
       -sforbrr(lsbeam)-sfdtbmc(lsbeam)
  !
  !  cx fast neutral source components
  snbxtot(lsbeam)=snbxtot(lsbeam)*zdti
  snbx1d(:,lsbeam)=snbx1d(:,lsbeam)*zdti
  snbxw0(lsbeam)=snbxw0(lsbeam)*zdti
  snbxv0(lsbeam)=snbxv0(lsbeam)*zdti
  snbxbb0(lsbeam)=snbxbb0(lsbeam)*zdti
  snbxbb1(lsbeam)=snbxbb1(lsbeam)*zdti

  !  cx recapture components
  sfxrcsc(lsbeam)=sfxrcsc(lsbeam)*zdti
  sfxrciz(lsbeam)=sfxrciz(lsbeam)*zdti
  sfxrccx(lsbeam)=sfxrccx(lsbeam)*zdti
  sfrbbiz(lsbeam)=sfrbbiz(lsbeam)*zdti
  sfrbbcx(lsbeam)=sfrbbcx(lsbeam)*zdti
  !
  return
  !
contains
  subroutine zwrite(zlbl,zoval,znval)
    character(len=*), intent(in) :: zlbl
    real(fp), intent(in) :: zoval,znval
    real(fp) :: zdif
    !
    zdif = abs(znval-zoval)/max(1.0e-30_fp,abs(znval),abs(zoval))
    write(nonlin,1001) zlbl, zoval, znval, zdif
1001 format(3x,a,' old=',1pe11.4,' new=',1pe11.4,' difn=',1pe11.4)
    !
    return
    !
  end subroutine zwrite
end subroutine outptb
!******************** end file outptb.for ; group outptb ***************
! 11sep2004 fgtok
