 Subroutine calc_sed_bed

    implicit none
    !
    !  Local variable declarations.
    !
    integer :: Ksed, i, ised,ibed, j, k, ks
    integer :: bnew, nnn, cnt,cc,jj

    real(rkind), parameter :: eps = 1.0E-14_rkind

    real(rkind) :: cff, cff1, cff2, cff3
    real(rkind) :: thck_avail, thck_to_add

    real(rkind), dimension(NST) :: dep_mass
    real(rkind), dimension(0:npa) :: tau_w

    real(rkind),allocatable,dimension(:,:)::temp

    if(myrank==0) write(*,*) "Start:  calc_sed_bed"


    if(BEDLOAD)then
      bnew=nnew
    else
      bnew=nstp
    endif

    !
    !-----------------------------------------------------------------------
    ! Compute sediment bed layer stratigraphy.
    !-----------------------------------------------------------------------
    !
    if(BEDLOAD_MPM.or.SUSLOAD)then
       tau_w = taub
    endif
    !
    !-----------------------------------------------------------------------
    !  Update bed properties according to ero_flux and dep_flux.
    !-----------------------------------------------------------------------
    !
    if(SUSLOAD)then
       NODE_LOOP : DO i=1,np
          SED_LOOP: DO ised=1,NST
             !
             !  The deposition and resurkindension of sediment on the bottom "bed"
             !  is due to precepitation flux FC(:,0), already computed, and the
             !  resurkindension (erosion, hence called ero_flux). The resurkindension is
             !  applied to the bottom-most grid box value qc(:,1) so the total mass
             !  is conserved. Restrict "ero_flux" so that "bed" cannot go negative
             !  after both fluxes are applied.
             !
             dep_mass(ised)=0.0_rkind

             if( SED_MORPH)then
                !
                ! Apply morphology factor.
                !
                sedbed%ero_flux(i,ised)=sedbed%ero_flux(i,ised)*sed(ised)%morph_fac
                sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)*            &
                     &                              sed(ised)%morph_fac
                !
             endif
             IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt.        &
                  &           0.0_rkind) THEN
                !
                !  If first time step of deposit, then store deposit material in
                !  temporary array, dep_mass.
                !
                IF ((T_model.gt.(bed(i,1,iaged)+1.1_rkind*DTsed)).and.   &
                     &            (bed(i,1,ithck).gt.newlayer_thick)) THEN
                   dep_mass(ised)=sedbed%settling_flux(i,ised)-               &
                        &                           sedbed%ero_flux(i,ised)
                END IF
                bed(i,1,iaged)=T_model
             END IF
             !
             !  Update bed mass arrays.
             !
             sedbed%bed_mass(i,1,nnew,ised)=MAX(sedbed%bed_mass(i,1,bnew,ised)-    &
                  &  (sedbed%ero_flux(i,ised)-           &
                  &  sedbed%settling_flux(i,ised)),     &
                  &  0.0_rkind)
             DO k=2,Nbed
                sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k,nstp,ised)
             END DO
          END DO SED_LOOP
          !
          !  If first time step of deposit, create new layer and combine bottom
          !  two bed layers.
          !
          cff=0.0_rkind
          !
          !  Determine if deposition ocurred here.
          !
          IF (Nbed.gt.1) THEN
             DO ised=1,NST
                cff=cff+dep_mass(ised)
             END DO
             IF (cff.gt.0.0_rkind) THEN
                !
                !  Combine bottom layers.
                !
                bed(i,Nbed,iporo)=0.5_rkind*(bed(i,Nbed-1,iporo)+        &
                     &                                    bed(i,Nbed,iporo))
                bed(i,Nbed,iaged)=0.5_rkind*(bed(i,Nbed-1,iaged)+        &
                     &                                    bed(i,Nbed,iaged))
                DO ised=1,NST
                   sedbed%bed_mass(i,Nbed,nnew,ised)=                           &
                        &  sedbed%bed_mass(i,Nbed-1,nnew,ised)+      &
                        &  sedbed%bed_mass(i,Nbed  ,nnew,ised)
                END DO
                !
                !  Push layers down.
                !
                DO k=Nbed-1,2,-1
                   bed(i,k,iporo)=bed(i,k-1,iporo)
                   bed(i,k,iaged)=bed(i,k-1,iaged)
                   DO ised =1,NST
                      sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k-1,nnew,ised)
                   END DO
                END DO
                !
                !  Set new top layer parameters.
                !
                DO ised=1,NST
                   sedbed%bed_mass(i,2,nnew,ised)=MAX(sedbed%bed_mass(i,2,nnew,ised)-&
                        &  dep_mass(ised),0.0_rkind)
                   sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised)
                END DO
             END IF
          END IF !NBED=1
          !
          ! Recalculate thickness and fractions for all layers.
          !
          DO k=1,Nbed
             cff3=0.0_rkind
             DO ised=1,NST
                cff3=cff3+sedbed%bed_mass(i,k,nnew,ised)
             END DO
             IF (cff3.eq.0.0_rkind) THEN
                cff3=eps
             END IF
             bed(i,k,ithck)=0.0_rkind
             DO ised=1,NST
                sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3
                bed(i,k,ithck)=MAX(bed(i,k,ithck)+                    &
                     &  sedbed%bed_mass(i,k,nnew,ised)/               &
                     &  (sed(ised)%Srho*                          &
                     &  (1.0_rkind-bed(i,k,iporo))),0.0_rkind)
             END DO
          END DO
       END DO NODE_LOOP
       !
       !  End of suspended Sediment only section.
       !
    endif
    !
    !  Ensure top bed layer thickness is greater or equal than active layer
    !  thickness. If need to add sed to top layer, then entrain from lower
    !  levels. Create new layers at bottom to maintain Nbed.
    !
    NODE_LOOP2 : DO i=1,np
       !
       !  Calculate active layer thickness, bottom(i,iactv).
       !
       bottom(i,iactv)=MAX(0.0_rkind,                                 &
            &   0.007_rkind*                               &
            &   (tau_w(i)-bottom(i,itauc)))+   &
            &    6.0_rkind*bottom(i,isd50)
       !
       if (SED_MORPH)then
          !
          ! Apply morphology factor.
          !
          bottom(i,iactv)=MAX(bottom(i,iactv)*sed(1)%morph_fac,      &
               &   bottom(i,iactv))
       endif
       !
       IF (bottom(i,iactv).gt.bed(i,1,ithck)) THEN
          IF (Nbed.eq.1) THEN
             bottom(i,iactv)=bed(i,1,ithck)
          ELSE
             thck_to_add=bottom(i,iactv)-bed(i,1,ithck)
             thck_avail=0.0_rkind
             Ksed=1                                        ! initialize
             DO k=2,Nbed
                IF (thck_avail.lt.thck_to_add) THEN
                   thck_avail=thck_avail+bed(i,k,ithck)
                   Ksed=k
                END IF
             END DO
             !
             !  Catch here if there was not enough bed material.
             !
             IF (thck_avail.lt.thck_to_add) THEN
                bottom(i,iactv)=bed(i,1,ithck)+thck_avail
                thck_to_add=thck_avail
             END IF
             !
             !  Update bed mass of top layer and fractional layer.
             !
             cff2=MAX(thck_avail-thck_to_add,0.0_rkind)/                  &
                  &   MAX(bed(i,Ksed,ithck),eps)
             DO ised=1,NST
                cff1=0.0_rkind
                DO k=1,Ksed
                   cff1=cff1+sedbed%bed_mass(i,k,nnew,ised)
                END DO
                cff3=cff2*sedbed%bed_mass(i,Ksed,nnew,ised)
                sedbed%bed_mass(i,1   ,nnew,ised)=cff1-cff3
                sedbed%bed_mass(i,Ksed,nnew,ised)=cff3
             END DO
             !
             !  Update thickness of fractional layer ksource_sed.
             !
             bed(i,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0_rkind)
             !
             !  Upate bed fraction of top layer.
             !
             cff3=0.0_rkind
             DO ised=1,NST
                cff3=cff3+sedbed%bed_mass(i,1,nnew,ised)
             END DO
             IF (cff3.eq.0.0_rkind) THEN
                cff3=eps
             END IF
             DO ised=1,NST
                sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3
             END DO
             !
             !  Upate bed thickness of top layer.
             !
             bed(i,1,ithck)=bottom(i,iactv)
             !
             !  Pull all layers closer to the surface.
             !
             DO k=Ksed,Nbed
                ks=Ksed-2
                bed(i,k-ks,ithck)=bed(i,k,ithck)
                bed(i,k-ks,iporo)=bed(i,k,iporo)
                bed(i,k-ks,iaged)=bed(i,k,iaged)
                DO ised=1,NST
                   sedbed%bed_frac(i,k-ks,ised)=sedbed%bed_frac(i,k,ised)
                   sedbed%bed_mass(i,k-ks,nnew,ised)=sedbed%bed_mass(i,k,nnew,ised)
                END DO
             END DO
             !
             !  Add new layers onto the bottom. rkindlit what was in the bottom layer to
             !  fill these new empty cells. ("ks" is the number of new layers).
             !
             ks=Ksed-2
             cff=1.0_rkind/REAL(ks+1,rkind)
             DO k=Nbed,Nbed-ks,-1
                bed(i,k,ithck)=bed(i,Nbed-ks,ithck)*cff
                bed(i,k,iaged)=bed(i,Nbed-ks,iaged)
                DO ised=1,NST
                   sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised)
                   sedbed%bed_mass(i,k,nnew,ised)=                            &
                        &  sedbed%bed_mass(i,Nbed-ks,nnew,ised)*cff
                END DO
             END DO
          END IF  ! Nbed > 1
       END IF  ! increase top bed layer

    END DO NODE_LOOP2
    !
    !-----------------------------------------------------------------------
    ! Store old bed thickness.
    !-----------------------------------------------------------------------
    !
    if (SED_MORPH)then
       do i=1,np
          sedbed%bed_thick(i,nnew)=0.0_rkind
          DO k=1,Nbed
             sedbed%bed_thick(i,nnew)=sedbed%bed_thick(i,nnew)+                  &
                  &                       bed(i,k,ithck)
          END DO
       end do

       call exchange_p2d(sedbed%bed_thick(:,nnew))
 
       bottom(i,morph) = sedbed%bed_thick(i,nnew)
    endif


       allocate(temp(0:npa,1:nbed))
       DO ised=1,NST
          temp = 0.0_rkind
          temp = sedbed%bed_frac(:,:,ised)
          !call aexchange(nc,myid,nprocs,temp)
          sedbed%bed_frac(:,:,ised) = temp
          temp = 0.0_rkind
          temp = sedbed%bed_mass(:,:,nnew,ised)
          !call aexchange(nc,myid,nprocs,temp)
          sedbed%bed_mass(:,:,nnew,ised) = temp
       end do

       do ibed=1,MBEDP
          temp = 0.0_rkind      
          temp = bed(:,:,ibed)
         ! call aexchange(nc,myid,nprocs,temp)
          bed(:,:,ibed)=temp
       end do
       deallocate(temp)

    if(myrank==0) write(*,*) "End:  calc_sed_bed"

    return   

  End Subroutine calc_sed_bed