!******************** start file depnrm.for ; group depnrm *************
!---------------------------------------------------------------------
!  depnrm - renormalize dep output profiles - cted from old nbmain
!   dmc mar 1985
!   passed quantities:
!     zwgrat - power renormalizer
!     zwrats - ptcl/ sce fcn renormalizer
!
!    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 depnrm should also be smoothed
!    in subroutine depsmoo **
!    see also:  outptb2, orbsmoo2.
!    see also:  outptb, orbsmoo.
!
!  generally the time normalization (1/(beam timestep))
!  is applied here to quantities which are strictly dep output
!  (bn00; dep only source fcns)
!
!  the time normalization is not applied to quantities which will
!  subsequently have contributions added due to nth generation cx
!  (i.e. nutrav)
!
!  ** see also outptb and outptb2 subroutines -- finally renormalization
!  of outputs of the beam monte carlo code
!
!  **bugfix dmc sept 2006**
!  when ndep0_use > ndep some quantities did not get the right normalization
!    (quantities that get sum contributions only in the lower part of dep
!    where new particles are created:  bphdeps & all nb0jxb output)
!    apply zwgfac to correct!
!
!-------------------------------------------
!  dmc aug 1988 -- multispecies beam code conversion
!
!   depnrm now gets called n times, once for each specie of beam
!  being deposited.  common index lsbeam indicates the current
!  beam specie.
!
!   scalar and 1d profile quantities are now generally stored in
!  two ways-- one array indexed by beam species and another array
!  or scalar giving the sum over all species.
!
!   2d and fancier profiles are generally stored indexed by beam
!  species but the sum over all species is not provided.
!
!   depnrm will compute the required sums
!
!  **depnrm gets called for all beam species** even if no new
!  ptcls are deposited for a given species.
!
!-------------------------------------------
!
subroutine depnrm(zwgrat,zwrats,zwrats_tr,zwgfac)
  !
  use nbi_com
  use nbi_mpi_com
  !
  implicit none
  real(fp) :: bsumint,zwnorm,zphjxb
  !
  real(fp) :: zn00t(3,mj)
  !
  !  zwgrat -- energy based renormalization
  !  zwrats -- particle based renormalization
  !  zwgfac -- enhancement for sums not incremented in ndep0_use loop
  !
  real(fp) :: zwgrat,zwrats,zwgfac,zwrats_tr
  !
  integer :: j,jp1,ib,ien,ig,ifun,ifi,ifn,isb,iz,ino,i,idir,ia,ie
  integer :: ish,jiav,ix,inum
  integer :: iindf,iioffn,iindn,jk
  !
  real(fp) :: zdti,znorm,zvnrm,zvnrm0,zchk,znchk,zphchk,zbpl,zpchk
  real(fp) :: zn0mx  ! (avoid compiler warning)
  logical :: key_myidd_0
  !
  !------------------------------------------------------------------
  !
  !  ztdep common replaced by nbi_com variables
  !
  !  common block also used in dep -- testing fpp first orbit loss
  !    accounting
  !   xn_dep,xn_orb-- #/sec deposited, left orbitting
  !   xp_dep,xp_orb-- power deposited, left orbitting
  !   xph_dep,xph_orb-- ang. mo./sec deposited, left orbitting
  !
  !   differences to be accounted by:  prompt orbit losses
  !   jxb torque and associated work term
  !
  !------------------------------------------------------------------
  !
  zdti=1.0_fp/delt
  !
  !  survival fraction array...
  !
  do ib=1,mib
     do ien=1,3
        inum=wfrac_survive(ien,ib)+0.5_fp
        if(inum.eq.0) then
           frac_survive(ien,ib)=1.0_fp
        else
           frac_survive(ien,ib)=sfrac_survive(ien,ib)/inum
        end if
     end do
  end do
  !
  key_myidd_0=.false.
  if(myidd.eq.0) key_myidd_0=.true.
#ifdef __MPI
#else
  !to allow  serial code when debugging mpi
  key_myidd_0=.true.
#endif
  if(key_myidd_0) then
     if(nltrk_dep0.and.mibs2use(lsbeam).gt.0) then
        track_sinj(1:ndep0_use,mibs2use(lsbeam)) = zwrats_tr* &
             track_sinj(1:ndep0_use,mibs2use(lsbeam))
     end if
  end if
  !
  !  scalars:
  !
  !     russian-roulette -- real weight of old ions, do not use zwrats!
  !
  sfdeprr(lsbeam)=sfdeprr(lsbeam)*zdti ! ions displaced for new deposition
  !
  sforbrr(lsbeam)=sforbrr(lsbeam)*zdti ! underweight ions
  !
  !   shinethru power
  !
  bpshins(lsbeam)=bpshins(lsbeam)*zwgrat*xe*zdti
  !
  sfshine(lsbeam)=sfshine(lsbeam)*zwrats*zdti
  !
  !   integrated deposition sources, by category
  sfdepsc(lsbeam)=sfdepsc(lsbeam)*zwrats*zdti  ! total
  sfdepiz(lsbeam)=sfdepiz(lsbeam)*zwrats*zdti  ! th. ionization
  sfdepcx(lsbeam)=sfdepcx(lsbeam)*zwrats*zdti  ! th. cx
  sfdbbcx(lsbeam)=sfdbbcx(lsbeam)*zwrats*zdti  ! beam-beam cx
  sfdbbiz(lsbeam)=sfdbbiz(lsbeam)*zwrats*zdti  ! beam-beam ioniz.

  !
  !  jxb:  flr->g.c.; may also include 1st orbit effect if orbits
  !  done by fpp
  !
  bpjxbs_dep(lsbeam)=bpjxbs_dep(lsbeam)*zwrats*zwgfac  ! dmc sept 2006
  bpjxbs(lsbeam)=0.0_fp             ! sum reused, nb0jxb in orbit loop
  !
  !  to limiter (if dep doing first orbit calculations)
  !
  bpls_dep(lsbeam)=bpls_dep(lsbeam)*zwgrat        ! power
  !c
  bphorbs_dep(lsbeam)=bphorbs_dep(lsbeam)*zwrats  ! momentum
  !
  bxflws_dep(lsbeam)=bxflws_dep(lsbeam)*zwrats    ! particles

  !
  !   angular momentum deposition
  !
  bphdeps(lsbeam)=bphdeps(lsbeam)*zwrats*zwgfac   ! dmc sept 2006, at point of ionization
  bphdepsgc(lsbeam)=bphdepsgc(lsbeam)*zwrats*zwgfac   !  at gc
  bphinjs(lsbeam)=bphinjs(lsbeam)*zwrats   ! injected momentum
  bphdepsgc(lsbeam)=bphdepsgc(lsbeam)*zwrats*zwgfac   !  at gc
  bphinjs(lsbeam)=bphinjs(lsbeam)*zwrats   ! injected momentum
  !
  !     shine through momentum
  bphshins(lsbeam)=bphshins(lsbeam)*zwrats
  !
  !  jxb:  flr->g.c.
  bphjxbs_dep(lsbeam)=bphjxbs_dep(lsbeam)*zwrats*zwgfac  ! dmc sept 2006
  bphjxbs(lsbeam)=0.0_fp            ! sum reused, nb0jxb in orbit loop
  !---------------------------------------
  !  shinethru and prompt orbit losses, by beam sce index
  !  only renormalize sces of the current species
  !
  if(.not.nlfprod(lsbeam)) then
     do ib=1,nbeam
        if(lsbeam.eq.lsbeama(ib)) then
           bpshna(ib,1:3)=bpshna(ib,1:3)*zwgrat*xe*zdti
           bpprma(ib)=bpprma(ib)*zwgrat*xe*zdti
        end if
     end do
  end if

  iioffn=6
  iindf=1
  iindn=iindf+iioffn

  do jk = lcentr,ledge
     taubsl(jk,iindf,lsbeam)=taubsl(jk,iindf,lsbeam) &
          +taubsl_tmp(jk,iindf,lsbeam)
     taubsl(jk,iindn,lsbeam)=taubsl(jk,iindn,lsbeam) &
          +taubsl_tmp(jk,iindn,lsbeam)

     taubpa(jk,iindf,lsbeam)=taubpa(jk,iindf,lsbeam) &
          +taubpa_tmp(jk,iindf,lsbeam)
     taubpa(jk,iindn,lsbeam)=taubpa(jk,iindn,lsbeam) &
          +taubpa_tmp(jk,iindn,lsbeam)

     mc_newdens(jk,lsbeam) = mc_newdens(jk,lsbeam)*zwrats/dvolzc(jk) &
          *zwgfac
  end do
  !
  !---------------------------------------
  !
  !  1d profile outputs
  !
  !  includes deposition-only plotting outputs
  !  other terms will have nth generation cx contributions added in
  !
  !  as with all beam code outputs, there is *no* volume normalization
  !
  !  also estimate beam-beam ch-ex and impact ionz. sce/sink fcns
  !  from previous timestep arrays
  !

  !      replaced resetr call with f90 assign jed 8/01
  !      call resetr(zn00t,3*mj,0.0)

  zn00t=0.0_fp

  !dbg
  !dbg	call zsdepchk(-1)
  !dbg
  do j=lcentr,ledge
     jp1=j+1
     !
     sorbes_dep(j,lsbeam)=sorbes_dep(j,lsbeam)*zwrats
     !
     !  normalize 1st generation beam neutral density array
     !  calculated in dep; indexbed by energy (full,half,third)
     !
     sdep_oa(j,lsbeam)=sdep_oa(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=0.0_fp
     !
     sbbidps(j,lsbeam)=sbbidps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbbidps(j,lsbeam)
     !
     sbbxdps(j,lsbeam)=sbbxdps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbbxdps(j,lsbeam)
     !
     sbedeps(j,lsbeam)=sbedeps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbedeps(j,lsbeam)
     !
     sbiedeps(j,lsbeam)=sbiedeps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbiedeps(j,lsbeam)
     !
     sbizdeps(j,lsbeam)=sbizdeps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbizdeps(j,lsbeam)
     !
     sbildeps(j,lsbeam)=sbildeps(j,lsbeam)*zwrats*zdti
     bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbildeps(j,lsbeam)
     !
     do ien=1,3
        bdepe(j,ien,lsbeam)=0.0_fp
        do ig=1,ng
           sbcxdpe(j,ig,ien,lsbeam)=sbcxdpe(j,ig,ien,lsbeam)* &
                zwrats*zdti
           bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
                sbcxdpe(j,ig,ien,lsbeam)
        end do
        !
        sbbidpe(j,ien,lsbeam)=sbbidpe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbbidpe(j,ien,lsbeam)
        !
        sbbxdpe(j,ien,lsbeam)=sbbxdpe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbbxdpe(j,ien,lsbeam)
        !
        sbedepe(j,ien,lsbeam)=sbedepe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbedepe(j,ien,lsbeam)
        !
        sbiedepe(j,ien,lsbeam)=sbiedepe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbiedepe(j,ien,lsbeam)
        !
        sbizdepe(j,ien,lsbeam)=sbizdepe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbizdepe(j,ien,lsbeam)
        !
        sbildepe(j,ien,lsbeam)=sbildepe(j,ien,lsbeam)*zwrats*zdti
        bdepe(j,ien,lsbeam)=bdepe(j,ien,lsbeam)+ &
             sbildepe(j,ien,lsbeam)
     end do
     !
     if(.not.nlfprod(lsbeam)) then
        do ib=1,nbeam
           if(lsbeama(ib).eq.lsbeam) then
              do ien=1,3
                 if(bn00(ib,ien,j).gt.0.0_fp) vpvn00(ib,ien,j)= &
                      vpvn00(ib,ien,j)/bn00(ib,ien,j)
                 bn00(ib,ien,j)=bn00(ib,ien,j)*zwrats*zdti
                 bn00t(ien,j)=bn00t(ien,j)+bn00(ib,ien,j)
                 zn00t(ien,j)=zn00t(ien,j)+bn00(ib,ien,j)
                 bdepb(j,ien,ib)=bdepb(j,ien,ib)*zwrats*zdti

              end do
              !mg
              tqjxb_nb_dep(j,ib)=tqjxb_nb_dep(j,ib)*zwrats*zwgfac
              tqjxb_nb(j,ib)=0.0_fp ! sum reused, nb0jxb in orbit loop

              !mg

           end if

        end do
     end if
     !
     do ig=1,ng
        sorbhs_dep(j,ig,lsbeam)=sorbhs_dep(j,ig,lsbeam)*zwrats
        sorbn0s_dep(j,ig,lsbeam)=sorbn0s_dep(j,ig,lsbeam)*zwrats
        sbcxdps(j,ig,lsbeam)=sorbn0s_dep(j,ig,lsbeam)*zdti
        bdeps(j,lsbeam)=bdeps(j,lsbeam)+sbcxdps(j,ig,lsbeam)

     end do
     !
     tqbcxs_dep(j,lsbeam)=tqbcxs_dep(j,lsbeam)*zwrats
     !
     tqjxbs_dep(j,lsbeam)=tqjxbs_dep(j,lsbeam)*zwrats*zwgfac ! dmc sept 2006
     tqjxbs(j,lsbeam)=0.0_fp ! sum reused, nb0jxb in orbit loop
     !


     pbephis_dep(j,lsbeam)=pbephis_dep(j,lsbeam)*zwrats*zwgfac ! dmc sept 2006
     pbephis(j,lsbeam)=0.0_fp ! sum reused, nb0jxb in orbit loop
     !
     ! radial current contributions
     !
     tirtots_dep(jp1,lsbeam)=tirtots_dep(jp1,lsbeam)*zwrats*zwgfac ! dmc sept 06
     tirtots(jp1,lsbeam)=0.0_fp ! sum reused, nb0jxb in orbit loop
     fpixbts_dep(jp1,lsbeam)=fpixbts_dep(jp1,lsbeam)*zwrats*zwgfac ! dmc sept 06
     fpixbts(jp1,lsbeam)=0.0_fp ! sum reused, nb0jxb in orbit loop
     !
     ! dep tau output normalization
     taubsl(j,19,lsbeam)=0.0_fp
     taubsl(j,20,lsbeam)=0.0_fp
     taubpa(j,19,lsbeam)=0.0_fp
     taubpa(j,20,lsbeam)=0.0_fp
     ifun=6
     do ifi=1,ifun
        ifn=ifi+6
        taubsl(j,19,lsbeam)=taubsl(j,19,lsbeam)+taubsl(j,ifi,lsbeam)
        taubsl(j,20,lsbeam)=taubsl(j,20,lsbeam)+taubsl(j,ifn,lsbeam)
        if(taubsl(j,ifn,lsbeam).eq.0.0_fp) then
           taubsl(j,ifi,lsbeam)=0.0_fp
        else
           taubsl(j,ifi,lsbeam)= &
                taubsl(j,ifi,lsbeam)/taubsl(j,ifn,lsbeam)
        end if
        taubpa(j,19,lsbeam)=taubpa(j,19,lsbeam)+taubpa(j,ifi,lsbeam)
        taubpa(j,20,lsbeam)=taubpa(j,20,lsbeam)+taubpa(j,ifn,lsbeam)
        if(taubpa(j,ifn,lsbeam).eq.0.0_fp) then
           taubpa(j,ifi,lsbeam)=0.0_fp
        else
           taubpa(j,ifi,lsbeam)= &
                taubpa(j,ifi,lsbeam)/taubpa(j,ifn,lsbeam)
        end if

     end do
     if(taubsl(j,20,lsbeam).eq.0.0_fp) then
        taubsl(j,19,lsbeam)=0.0_fp
     else
        taubsl(j,19,lsbeam)=taubsl(j,19,lsbeam)/taubsl(j,20,lsbeam)
     end if
     if(taubpa(j,20,lsbeam).eq.0.0_fp) then
        taubpa(j,19,lsbeam)=0.0_fp
     else
        taubpa(j,19,lsbeam)=taubpa(j,19,lsbeam)/taubpa(j,20,lsbeam)
     end if


  end do
  !
  !  excited states stuff ...
  !
  if(nsigexc_use.gt.0) then
     do ien=1,3
        do j=lcentr,ledge
           if(zn00t(ien,j).gt.0.0_fp) then
              znorm=zwrats*zdti/zn00t(ien,j)
              excit0(j,ien,lsbeam)=excit0(j,ien,lsbeam)*znorm
           end if
        end do
     end do
  else
     excit0(lcentr:ledge,1:3,lsbeam)=1.0_fp
  end if
  !
  !-----------------------------
  !  spatial 2d profile data for new code
  !   loop over beam dist fcn zones
  !
  if(nlfprod(lsbeam)) goto 205
  do iz=1,nfbznsi
     !
     !  1st gen beam neutral population
     do ino=1,2
        do ien=1,3
           do ib=1,nbeam
              if(lsbeama(ib).eq.lsbeam) then
                 if((nsigexc_use.gt.0).and. &
                      (bn002(ib,ien,ino,iz).gt.0.0_fp)) then
                    do isb=1,nsbeam
                       hexbbcx0(ib,ien,ino,iz,isb)= &
                            hexbbcx0(ib,ien,ino,iz,isb)/ &
                            bn002(ib,ien,ino,iz)
                    end do
                 else
                    do isb=1,nsbeam
                       hexbbcx0(ib,ien,ino,iz,isb)=1.0_fp
                    end do
                 end if
                 bn002(ib,ien,ino,iz)=bn002(ib,ien,ino,iz)*zwrats* &
                      zdti
              end if         ! lsbeama test
           end do            ! ib
        end do               ! ien
        !
        !  renormalize monte carlo neutral velocity direction cosines sum
        !  (by sum**2=1)
        !
        do ib=1,nbeam
           if(lsbeama(ib).eq.lsbeam) then
              zvnrm=0.0_fp
              zn0mx=0.0_fp
              !     bugfix dmc - avoid floating underflow to zero from squaring small no.s
              zvnrm0=max( abs(vcxbn0(ib,1,ino,iz)), &
                   abs(vcxbn0(ib,2,ino,iz)), &
                   abs(vcxbn0(ib,3,ino,iz)))
              !
              !     check that there is non-zero neutral density
              do i=1,3
                 !     here i is energy index
                 zn0mx=max(zn0mx,bn002(ib,i,ino,iz))
                 !     here i is direction index
                 if(zvnrm0.gt.0.0_fp) zvnrm=zvnrm+ &
                      (vcxbn0(ib,i,ino,iz)/zvnrm0)**2
              end do
              !
              if(zn0mx.gt.0.0_fp) then
                 zvnrm=sqrt(zvnrm)*zvnrm0
                 !     neutral velocity avg direction cosine (3 directions idir)
                 do idir=1,3
                    vcxbn0(ib,idir,ino,iz)=vcxbn0(ib,idir,ino,iz) &
                         /zvnrm
                 end do
              end if
           end if
           !
        end do               ! ib
        !
     end do                  ! ino
  end do                     ! iz
  !
  !     thermal neutral source fcns
  do ig=1,ng
     do iz=1,nfbzns
        srbn02_dep(iz,ig,lsbeam)=srbn02_dep(iz,ig,lsbeam)*zwrats
     end do
  end do
  !
205 continue
  !--------------------------------------------------------------------
  !  for the fpp code  -- debug code also
  !
  if(nlfbmfpp(lsbeam)) then
     !
     zchk=0.0_fp
     do jiav=lcentr,lep1
        fbdepco(jiav,lsbeam)=fbdepco(jiav,lsbeam)*zwgrat*zdti
        fbdepct(jiav,lsbeam)=fbdepct(jiav,lsbeam)*zwgrat*zdti
     end do
  end if
  !
  !  internal ptcl balance check
  if(kdebug1.eq.1) then
     xn_dep(lsbeam)=xn_dep(lsbeam)*zwrats
     xn_orb(lsbeam)=xn_orb(lsbeam)*zwrats
     znchk=xn_dep(lsbeam)-xn_orb(lsbeam)-bxflws_dep(lsbeam)
     write(nonlin,3001) xn_dep(lsbeam),bxflws_dep(lsbeam), &
          xn_orb(lsbeam),zchk,znchk
3001 format(/' dep ptcl bal check:'/ &
          '  ionized:  ',1pe12.5,'  1st orbit loss:  ',1pe12.5/ &
          '  orbiting: ',1pe12.5,'  same in fb(fpp): ',1pe12.5/ &
          '  bal. chk: ',1pe12.5)
     !
     !     internal momentum balance check
     xph_dep(lsbeam)=xph_dep(lsbeam)*zwrats
     xph_orb(lsbeam)=xph_orb(lsbeam)*zwrats
     zphjxb=bsumint(tqjxbs_dep(1,lsbeam))
     zphchk=xph_dep(lsbeam)-zphjxb &
          -bphorbs_dep(lsbeam)-xph_orb(lsbeam)
     write(nonlin,3003) xph_dep(lsbeam),bphorbs_dep(lsbeam), &
          zphjxb,xph_orb(lsbeam),zphchk
3003 format(/' dep momentum bal check:'/ &
          '  ionized:  ',1pe12.5,'  1st orbit loss:  ',1pe12.5/ &
          '  jxb tq:   ',1pe12.5,'  left orbiting:   ',1pe12.5/ &
          '  bal. chk: ',1pe12.5)
     !
     !     internal energy balance check
     xp_dep(lsbeam)=xp_dep(lsbeam)*zwrats
     xp_orb(lsbeam)=xp_orb(lsbeam)*zwrats
     zbpl=bpls_dep(lsbeam)*(zwrats/zwgrat)
     zpchk=xp_dep(lsbeam)-bpjxbs_dep(lsbeam) &
          -zbpl-xp_orb(lsbeam)
     write(nonlin,3005) xp_dep(lsbeam),zbpl, &
          bpjxbs_dep(lsbeam),xp_orb(lsbeam),zpchk
3005 format(/' dep power bal check:'/ &
          '  ionized:  ',1pe12.5,'  1st orbit loss:  ',1pe12.5/ &
          '  jxb tq:   ',1pe12.5,'  left orbiting:   ',1pe12.5/ &
          '  bal. chk: ',1pe12.5)
     !
     write(nonlin,3007) zwgrat,zwrats
3007 format(/ &
          ' dep:  pwr renorm:  ',1pe12.5,'  sce renorm:',1pe12.5)
     !
  end if
  !
  !     renorm trapping fraction calculation
  !
  zwnorm=bsumint(wtrap0(1:mj,lsbeam))
  if(zwnorm.gt.0.0_fp) then
     tbtrap0(lsbeam)=bsumint(btrap0(1:mj,lsbeam))/zwnorm
  end if
  !
  do j=lcentr,lep1
     if(wtrap0(j,lsbeam).ne.0.0_fp) then
        btrap0(j,lsbeam)=btrap0(j,lsbeam)/wtrap0(j,lsbeam)
     else
        btrap0(j,lsbeam)=-0.1_fp
     end if
  end do
  !
  !     dmc mar. 2002:  re-organized tau(pa) and tau(sl) output profiles...
  !
  taubsl0(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,19,lsbeam)
  taubpa0(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,19,lsbeam)
  if(.not.nlfprod(lsbeam)) then
     ! co
     taubslco1(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,1,lsbeam)
     taubslco2(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,2,lsbeam)
     taubslco3(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,3,lsbeam)
     ! ctr
     taubslctr1(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,4,lsbeam)
     taubslctr2(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,5,lsbeam)
     taubslctr3(lcentr:ledge,lsbeam)=taubsl(lcentr:ledge,6,lsbeam)
     ! co
     taubpaco1(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,1,lsbeam)
     taubpaco2(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,2,lsbeam)
     taubpaco3(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,3,lsbeam)
     ! ctr
     taubpactr1(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,4,lsbeam)
     taubpactr2(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,5,lsbeam)
     taubpactr3(lcentr:ledge,lsbeam)=taubpa(lcentr:ledge,6,lsbeam)
  end if

  return
end subroutine depnrm
!******************** end file depnrm.for ; group depnrm ***************
