 Subroutine calc_sed_bed2

    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) :: nlysm

    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_bed2"

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

    ! KLUDGE alert
    ! minlayer_thick = 0.0005
      minlayer_thick = newlayer_thick
    !
    !-----------------------------------------------------------------------
    ! 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
          !
          !  The deposition and resurkindension of sediment on the bottom "bed"
          !  is due to precipitation 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.
          !
          SED_LOOP: DO ised=1,NST  
             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

             !  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

          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
          bed(i,1,ithck)=0.0_rkind
          DO ised=1,NST
             sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3
             bed(i,1,ithck)=MAX(bed(i,1,ithck)+                      &
                  &     sedbed%bed_mass(i,1,nnew,ised)/                                &
                  &     (sed(ised)%Srho*                                         &
                  &     (1.0_rkind-bed(i,1,iporo))),0.0_rkind)
          END DO

       END DO NODE_LOOP
    endif  !/* SUSPLOAD section */
    !
    !-----------------------------------------------------------------------
    !  At this point, all deposition or erosion is complete, and
    !  has been added/subtracted to top layer. Thickness has NOT been corrected.
    !-----------------------------------------------------------------------
    !
    NODE_LOOP2 : DO i=1,np
       !   Calculate active layer thickness, bottom(i,iactv).
       !   (trunk version allows this to be zero...this has minimum of 6*D50)
       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
       !
       !   Calculate net deposition and erosion
       cff=0.0_rkind
       cff2=0.0_rkind
       DO ised=1,NST
          cff=cff+sedbed%settling_flux(i,ised)
          cff2=cff2+sedbed%ero_flux(i,ised)
          dep_mass(ised)=0.0_rkind
          IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt.      &
               &     0.0_rkind) THEN
             dep_mass(ised)=sedbed%settling_flux(i,ised)-              &
                  &        sedbed%ero_flux(i,ised)
          END IF
       END DO

       IF ( cff-cff2.GT.0.0_rkind) THEN ! NET depostion
          !  Deposition. Determine if we need to create a new bed layer 
          ! (no test for age here)
          bed(i,1,iaged)=T_model
          IF(bed(i,1,ithck).gt.                                   &
               &    MAX(bottom(i,iactv),newlayer_thick)) THEN
             ! Top layer is too thick
             IF (Nbed.gt.2) THEN
                IF(bed(i,2,ithck).lt.minlayer_thick) THEN
                   ! Layer 2 is smaller than minimum size
                   ! Instead of pushing down all layers, just combine top 2 layers
                   cff=0.0_rkind
                   cff1=0.0_rkind
                   cff2=0.0_rkind
                   DO ised=1,NST
                      cff =cff +dep_mass(ised)
                      cff1=cff1+sedbed%bed_mass(i,1,nnew,ised)
                      cff2=cff2+sedbed%bed_mass(i,2,nnew,ised)
                   END DO
                   !  Update bed mass
                   DO ised=1,NST
                      sedbed%bed_mass(i,2,nnew,ised)=                    &
                           &                 MAX(sedbed%bed_mass(i,2,nnew,ised)+           &
                           &                 sedbed%bed_mass(i,1,nnew,ised)-               &
                           &                 dep_mass(ised),0.0_rkind)
                      sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised)
                   END DO
                   ! ALA - average T_model and porosity
                   ! ALA CHECK WITH CRS cff1 or cff1-cff for first layer
                   bed(i,2,iaged)=(bed(i,1,iaged)*cff1+         &
                        &               bed(i,2,iaged)*cff2)/(cff1+cff2)
                   bed(i,1,iaged)=T_model
                   bed(i,2,iporo)=(bed(i,1,iporo)*cff1+         &
                        &               bed(i,2,iporo)*cff2)/(cff1+cff2)
                   ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER
                   bed(i,1,iporo)=bed(i,1,iporo)
                ELSE
                   ! Layer 2 is > minlayer thick, need another layer
                   !  Combine bottom layers.
                   cff1=0.0_rkind
                   cff2=0.0_rkind
                   DO ised=1,NST
                      cff1=cff1+sedbed%bed_mass(i,Nbed-1,nnew,ised)
                      cff2=cff2+sedbed%bed_mass(i,Nbed,nnew,ised)
                   END DO
                   bed(i,Nbed,iporo)=                             &
                        &              (bed(i,Nbed-1,iporo)*cff1+                &
                        &              bed(i,Nbed,iporo)*cff2)/(cff1+cff2)
                   bed(i,Nbed,iaged)=                             &
                        &              (bed(i,Nbed-1,iaged)*cff1+                &
                        &               bed(i,Nbed,iaged)*cff2)/(cff1+cff2)
                   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 parameters for top 2 layers
                   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
             ELSEIF (Nbed.eq.2) THEN 
                ! NBED=2
                cff1=0.0_rkind
                cff2=0.0_rkind
                DO ised=1,NST
                   cff1=cff1+sedbed%bed_mass(i,1,nnew,ised)
                   cff2=cff2+sedbed%bed_mass(i,2,nnew,ised)
                END DO
                DO ised=1,NST
                   sedbed%bed_mass(i,2,nnew,ised)=                       &
                        &            MAX(sedbed%bed_mass(i,2,nnew,ised)+              &
                        &                sedbed%bed_mass(i,1,nnew,ised)-              &
                        &                dep_mass(ised),0.0_rkind)
                   sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised)
                END DO
                ! ALA - average T_model and porosity
                bed(i,2,iaged)=(bed(i,1,iaged)*cff1+            &
                     &              bed(i,2,iaged)*cff2)/(cff1+cff2)
                bed(i,1,iaged)=T_model
                bed(i,2,iporo)=(bed(i,1,iporo)*cff1+            &
                     &              bed(i,2,iporo)*cff2)/(cff1+cff2)
                ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER
                bed(i,1,iporo)=bed(i,1,iporo)                    
             ELSE
                ! NBED=1
             END IF
          ELSE
             ! Net deposition has occured, but no new bed layer was created 
          END IF
       ELSE
          ! Net erosion occurred
          bed(i,1,iaged)=T_model
          IF (Nbed.eq.2) THEN 
             ! NBED=2
             DO ised=1,NST
                sedbed%bed_mass(i,2,nnew,ised)=                       &
                     &           MAX(sedbed%bed_mass(i,2,nnew,ised)+              &
                     &               sedbed%bed_mass(i,1,nnew,ised)-                  &
                     &               dep_mass(ised),0.0_rkind)
                sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised)
             END DO
          ELSEIF (Nbed.eq.1) THEN
             ! ALF NO NEED TO DO ANYTHING
          ELSE
          END IF
       END IF

       ! 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_LOOP2


    NODE_LOOP3 : DO i=1,np
       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)
             !
             !  Update 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
             ! ALA CHECK WITH CRS about bed_frac 
             nnn=0
             DO ised=1,NST
                nlysm(ised)=newlayer_thick*REAL(ks+1,rkind)*          &
                     &  (sed(ised)%Srho*                                 &
                     &  (1.0_rkind-bed(i,Nbed-ks,iporo)))*               &
                     &  sedbed%bed_frac(i,Nbed-ks,ised)
                IF (ks.gt.0) THEN
                   IF (sedbed%bed_mass(i,Nbed-ks,nnew,ised).gt.         &
                        & nlysm(ised)) THEN
                      nnn=nnn+1
                      nlysm(ised)=                                 &
                           &  newlayer_thick*REAL(ks,rkind)*        &
                           & (sed(ised)%Srho*                         &
                           &  (1.0_rkind-bed(i,Nbed-ks,iporo)))*       &
                           &  sedbed%bed_frac(i,Nbed-ks,ised)
                   END IF
                END IF
             END DO
             IF (nnn.eq.NST) THEN
                bed(i,Nbed,ithck)=bed(i,Nbed-ks,ithck)-           &
                     &  newlayer_thick*REAL(ks,rkind)
                DO ised=1,NST
                   sedbed%bed_mass(i,Nbed,nnew,ised)=                      &
                        &  sedbed%bed_mass(i,Nbed-ks,nnew,ised)-nlysm(ised)
                END DO
                DO k=Nbed-1,Nbed-ks,-1
                   bed(i,k,ithck)=newlayer_thick
                   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)=                      &
                           &  nlysm(ised)/REAL(ks,rkind)
                   END DO
                END DO
             ELSE
                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
          END IF  ! Nbed > 1
       END IF  ! increase top bed layer
    END DO NODE_LOOP3
    !
    !-----------------------------------------------------------------------
    ! 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
    !
    !-----------------------------------------------------------------------
    !  Apply periodic or gradient boundary conditions to property arrays.
    !-----------------------------------------------------------------------
    !
       allocate(temp(0:npa,1:nbed))
       DO ised=1,NST
          temp = 0.0_rkind
          temp = sedbed%bed_frac(:,:,ised)
          sedbed%bed_frac(:,:,ised) = temp
          temp = 0.0_rkind
          temp = sedbed%bed_mass(:,:,nnew,ised)
          sedbed%bed_mass(:,:,nnew,ised) = temp
       end do

       do ibed=1,MBEDP
          temp = 0.0_rkind      
          temp = bed(:,:,ibed)
          bed(:,:,ibed)=temp
       end do
       deallocate(temp)


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

    return

  End Subroutine calc_sed_bed2