#include "cppdefs.h"

      MODULE sed_flocs_mod

#if defined NONLINEAR && defined SEDIMENT && defined SUSPLOAD && defined SED_FLOCS

!
!svn $Id: sed_flocs.F 429 2009-12-20 17:30:26Z arango $
!=======================================================================
!  Copyright (c) 2002-2017 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license           Hernan G. Arango   !
!    See License_ROMS.txt                   Alexander F. Shchepetkin   !
!==================================================== John C. Warner ===
!                                                                      !
!  This routine computes floc transformation.                          !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Verney, R., Lafite, R., Claude Brun-Cottan, J., & Le Hir, P. (2011).!
!    Behaviour of a floc population during a tidal cycle: laboratory   !
!    experiments and numerical modelling. Continental Shelf Research,  !
!    31(10), S64-S83.                                                  !
!=======================================================================
!
      implicit none

      PRIVATE
      PUBLIC  :: sed_flocmod

      CONTAINS
!
!***********************************************************************
      SUBROUTINE sed_flocmod (ng, tile)
!***********************************************************************
!
      USE mod_param
      USE mod_forces
      USE mod_grid
      USE mod_mixing
      USE mod_ocean
      USE mod_stepping
      USE mod_bbl
      USE mod_sedflocs
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, iNLM, 16)
# endif
      CALL sed_flocmod_tile (ng, tile,                                 &
     &                       LBi, UBi, LBj, UBj, N(ng), NT(ng),        &
     &                       IminS, ImaxS, JminS, JmaxS,               &
     &                       nstp(ng), nnew(ng),                       &
     &                       GRID(ng) % z_r,                           &
     &                       GRID(ng) % z_w,                           &
     &                       GRID(ng) % Hz,                            &
# ifdef BBL_MODEL
     &                       BBL(ng) % bustrcwmax,                     &
     &                       BBL(ng) % bvstrcwmax,                     &
     &                       FORCES(ng) % Pwave_bot,                   &
# endif
     &                       FORCES(ng) % bustr,                       &
     &                       FORCES(ng) % bvstr,                       &
     &                       MIXING(ng) % Akt,                         &
     &                       MIXING(ng) % Akv,                         &
     &                       MIXING(ng) % Lscale,                      &
     &                       MIXING(ng) % gls,                         &
     &                       MIXING(ng) % tke,                         &
     &                       OCEAN(ng) % t,                            &
     &                       SEDFLOCS(ng) % f_mass,                    &
     &                       SEDFLOCS(ng) % f_diam)
# ifdef PROFILE
      CALL wclock_off (ng, iNLM, 16)
# endif
      RETURN
      END SUBROUTINE sed_flocmod
!
!***********************************************************************
      SUBROUTINE sed_flocmod_tile (ng, tile,                            &
     &                             LBi, UBi, LBj, UBj, UBk, UBt,        &
     &                             IminS, ImaxS, JminS, JmaxS,          &
     &                             nstp, nnew, z_r, z_w, Hz,            &
# ifdef BBL_MODEL
     &                             bustrcwmax,                          &
     &                             bvstrcwmax,                          &
     &                             Pwave_bot,                           &
# endif
     &                             bustr,                               &
     &                             bvstr,                               &
     &                             Akt,Akv,Lscale,gls,tke,              &
     &                             t,                                   &
     &                             f_mass,f_diam)
!***********************************************************************
!
      USE mod_param
      USE mod_scalars
      USE mod_sediment
      USE mod_sedflocs
!
      implicit none 
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: nstp, nnew
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: z_r(LBi:,LBj:,:)
      real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
      real(r8), intent(in) :: Hz(LBi:,LBj:,:)
#  ifdef BBL_MODEL
      real(r8), intent(in) :: bustrcwmax(LBi:,LBj:)
      real(r8), intent(in) :: bvstrcwmax(LBi:,LBj:)
      real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
#  endif
      real(r8), intent(in) :: bustr(LBi:,LBj:)
      real(r8), intent(in) :: bvstr(LBi:,LBj:)
      real(r8), intent(in) :: Akt(LBi:,LBj:,0:,:)
      real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
      real(r8), intent(in) :: Lscale(LBi:,LBj:,0:)
      real(r8), intent(in) :: tke(LBi:,LBj:,0:,:)
      real(r8), intent(in) :: gls(LBi:,LBj:,0:,:)
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) 
# else
      real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
#  ifdef BBL_MODEL
      real(r8), intent(in) :: bustrcwmax(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstrcwmax(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:UBk,3)
      real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:UBk,3)
      real(r8), intent(in) :: Lscale(LBi:UBi,LBj:UBj,0:UBk,3)
      real(r8), intent(in) :: tke(LBi:UBi,LBj:UBj,0:UBk,3)
      real(r8), intent(in) :: gls(LBi:UBi,LBj:UBj,0:UBk,3)
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
# endif
      real(r8), intent(in) :: f_mass(0:NCS+1)
      real(r8), intent(in) :: f_diam(NCS)
!
!  Local variable declarations.
! 
      integer :: i, indx, ised, j, k, ks
!
!  Variable declarations for floc model
!
      integer  :: iv1
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
      real(r8) :: Gval,diss,mneg,dttemp,f_dt
      real(r8) :: dt1,f_csum,epsilon8
      real(r8) :: cvtotmud,tke_av, gls_av, exp1, exp2, exp3, ustr2,effecz
      real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: susmud
      real(r8), dimension(N(ng),IminS:ImaxS,JminS:JmaxS) :: f_davg
      real(r8), dimension(N(ng),IminS:ImaxS,JminS:JmaxS) :: f_d50
      real(r8), dimension(N(ng),IminS:ImaxS,JminS:JmaxS) :: f_d90
      real(r8), dimension(N(ng),IminS:ImaxS,JminS:JmaxS) :: f_d10
      real(r8),dimension(1:NCS)     :: cv_tmp,NNin,NNout
!  f_mneg_param : negative mass tolerated to avoid small sub time step (g/l)
      real(r8), parameter :: f_mneg_param=0.000_r8
#include "set_bounds.h"

      epsilon8=epsilon(1.0)
!     epsilon8=1.e-8
!
!--------------------------------------------------------------------------
! * Executable part
!
      J_LOOP : DO j=Jstr,Jend
!
!  Extract mud variables from tracer arrays, place them into
!  scratch arrays, and restrict their values to be positive definite.
      DO k=1,N(ng)
        DO i=Istr,Iend
          Hz_inv(i,k)=1.0_r8/Hz(i,j,k)
        END DO
      END DO
      DO ised=1,NCS
         indx = idsed(ised)
         DO k=1,N(ng)
            DO i=Istr,Iend
!               susmud(i,k,ised)=MAX(t(i,j,k,nstp,indx),0.0_r8)
                susmud(i,k,ised)=t(i,j,k,nnew,indx)*Hz_inv(i,k)
            ENDDO
	 ENDDO
      ENDDO
 
! min concentration below which flocculation processes are not calculated
!      f_clim=0.001_r8 
       exp1 = 3.0_r8+gls_p(ng)/gls_n(ng)
       exp2 = 1.5_r8+gls_m(ng)/gls_n(ng)
       exp3 = -1.0_r8/gls_n(ng)

       DO i=Istr,Iend
           DO k=1,N(ng)
 
              f_dt=dt(ng)
              dttemp=0.0_r8

	     ! concentration of all mud classes in one grid cell
              cvtotmud=0.0_r8
              DO ised=1,NCS
                 cv_tmp(ised)=susmud(i,k,ised) 
                 cvtotmud=cvtotmud+cv_tmp(ised)

                 NNin(ised)=cv_tmp(ised)/f_mass(ised)
              ENDDO

              DO iv1=1,NCS
                 IF (NNin(iv1).lt.0.0_r8) THEN
                  WRITE(*,*) '***************************************'
                  WRITE(*,*) 'CAUTION, negative mass at cell i,j,k :',  &
     &                          i,j,k
                  WRITE(*,*) '***************************************'        
                 ENDIF
              ENDDO

              IF (cvtotmud .gt. f_clim) THEN

# if defined FLOC_TURB_DISS && !defined FLOC_BBL_DISS
!
!ALA dissipation from turbulence clossure
!
		 IF (k.eq.1) THEN
                    tke_av = tke(i,j,k-1,nnew)
                    gls_av = gls(i,j,k-1,nnew)
		 ELSEIF (k.eq.N(ng)) THEN
                    tke_av = tke(i,j,k,nnew)
                    gls_av = gls(i,j,k,nnew)
		 ELSE 
                    tke_av = 0.5_r8*(tke(i,j,k-1,nnew)+tke(i,j,k,nnew))
                    gls_av = 0.5_r8*(gls(i,j,k-1,nnew)+gls(i,j,k,nnew))
                ENDIF
!               exp1 = 3.0_r8+gls_p(ng)/gls_n(ng)
!               exp2 = 1.5_r8+gls_m(ng)/gls_n(ng)
!               exp3 = -1.0_r8/gls_n(ng)
                diss = gls_cmu0(ng)**exp1*tke_av**exp2*gls_av**exp3
# elif defined FLOC_BBL_DISS && !defined FLOC_TURB_DISS
!
!ALA dissipation from wavecurrent bottom stress
! NOT READY FOR PRIME TIME
! NEEDS VERTICAL DISTRIBUTION
! As first cut, use turbulence closure
!
		 IF (k.eq.1) THEN
                    tke_av = tke(i,j,k-1,nnew)
                    gls_av = gls(i,j,k-1,nnew)
		 ELSEIF (k.eq.N(ng)) THEN
                    tke_av = tke(i,j,k,nnew)
                    gls_av = gls(i,j,k,nnew)
		 ELSE 
                    tke_av = 0.5_r8*tke(i,j,k-1,nnew)+                  &
     &                       0.5_r8*tke(i,j,k,nnew)
                    gls_av = 0.5_r8*gls(i,j,k-1,nnew)+                  &
     &                       0.5_r8*gls(i,j,k,nnew)
                ENDIF
                diss = gls_cmu0(ng)**exp1*tke_av**exp2*gls_av**exp3
!
!  MODIFY THE BOTTOM LAYER TO INCLUDE WAVECURRENT STRESS
!
		IF (k.eq.1) THEN       
#  ifdef BBL_MODEL
                 ustr2 =sqrt((bustrcwmax(i,j)**2.0_r8+                  &
     &                       bvstrcwmax(i,j)**2.0_r8 ))
                 effecz = (ustr2**0.5_r8)*Pwave_bot(i,j)*0.5_r8/pi
                 diss = MAX((ustr2**(1.5_r8))/(vonKar*effecz),diss)
#  else
                 ustr2 =sqrt((bustr(i,j)**2.0_r8+bvstr(i,j)**2.0_r8 ))
                 diss = MAX((ustr2**(1.5_r8))/(vonKar*                  &
     &                       (z_r(i,j,1)-zw(i,j,0))),diss)
#  endif
               ENDIF
# else 
                diss = epsilon8
                IF (l_testcase) THEN
!                  if (j.eq.1.and.i.eq.1.and.k.eq.1) then
!                     WRITE(*,*) 'VERNEY ET AL TESTCASE FOR FLOCS'
!                  endif
                ELSE    
                  WRITE(*,*) 'CAUTION :'
                  WRITE(*,*) 'CHOOSE A DISSIPATION MODEL FOR FLOCS'
                  WRITE(*,*) 'SIMULATION STOPPED'
                  STOP
                ENDIF                 
# endif		           
                CALL flocmod_comp_g(k,i,j,Gval,diss,ng) 

                 DO WHILE (dttemp .le. dt(ng))

                    CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)
                    CALL flocmod_mass_control(NNout,mneg,ng)
                    IF (mneg .gt. f_mneg_param) THEN
                       DO WHILE (mneg .gt. f_mneg_param)
                          f_dt=MIN(f_dt/2.0_r8,dt(ng)-dttemp)
                          IF (f_dt.lt.epsilon8) THEN
		             CALL flocmod_mass_redistribute(NNin,ng)
                             dttemp=dt(ng)
                             exit
                          ENDIF
                          CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)
                          CALL flocmod_mass_control(NNout,mneg,ng)   
                       ENDDO
                    ELSE

                       IF (f_dt.lt.dt(ng)) THEN
                          DO while (mneg .lt.f_mneg_param) 
                            IF (dttemp+f_dt .eq. dt(ng)) THEN
                           CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)
                              exit
                            ELSE
                              dt1=f_dt
                              f_dt=MIN(2.0_r8*f_dt,dt(ng)-dttemp)
                           CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)
                              CALL flocmod_mass_control(NNout,mneg,ng)
                              IF (mneg .gt. f_mneg_param) THEN 
                                f_dt=dt1
                           CALL flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)
                                exit
                              ENDIF
                            ENDIF
                          ENDDO
                       ENDIF
                    ENDIF
                    dttemp=dttemp+f_dt

                    NNin(:)=NNout(:) ! update new Floc size distribution
                    ! redistribute negative masses IF any on positive classes, 
                    ! depends on f_mneg_param
		    CALL flocmod_mass_redistribute(NNin,ng)

                    IF (dttemp.eq.dt(ng)) exit

                 ENDDO ! loop on full dt

              ENDIF ! only if cvtotmud > f_clim


              ! update mass concentration for all mud classes
              DO ised=1,NCS
                 susmud(i,k,ised)=NNin(ised)*f_mass(ised)
              ENDDO
           ENDDO
        ENDDO
!
!-----------------------------------------------------------------------
!  Update global tracer variables.
!-----------------------------------------------------------------------
!
      DO ised=1,NCS
         indx = idsed(ised)
          DO k=1,N(ng)
            DO i=Istr,Iend
              t(i,j,k,nnew,indx)=susmud(i,k,ised)*Hz(i,j,k)
            ENDDO
	 ENDDO
      ENDDO
      
      
      END DO J_LOOP
!      WRITE(*,*) 'END flocmod_main'
      END SUBROUTINE sed_flocmod_tile



!!===========================================================================
      SUBROUTINE flocmod_collfrag(Gval,ng,f_g4,f_l4)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_collfrag  ***
  !&E
  !&E ** Purpose : computation of collision fragmentation term, 
  !&E              based on McAnally and Mehta, 2001
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_comp_fsd
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
      USE mod_sedflocs
      USE mod_param
      USE mod_scalars
!
      implicit none 
      integer, intent(in) :: ng
      real(r8),intent(in) :: Gval

  !! * Local declarations
      integer      :: iv1,iv2,iv3
      real(r8) :: f_fp,f_fy,f_cfcst,gcolfragmin,gcolfragmax
      real(r8) :: gcolfragiv1,gcolfragiv2,f_weight,mult
      real(r8) :: cff1, cff2
      real(r8),DIMENSION(NCS,NCS,NCS) :: f_g4
      real(r8),DIMENSION(NCS,NCS) :: f_l4
      real(r8),DIMENSION(NCS) :: fdiam,fmass
! shorten names
      fdiam=SEDFLOCS(ng)%f_diam
      fmass=SEDFLOCS(ng)%f_mass
!  !!--------------------------------------------------------------------------
  !! * Executable part

      f_fp=0.1_r8
      f_fy=1e-10
      f_cfcst=3.0_r8/16.0_r8
      f_g4(1:NCS,1:NCS,1:NCS)=0.0_r8
      f_l4(1:NCS,1:NCS)=0.0_r8
      cff1=2.0_r8/(3.0_r8-f_nf)
      cff2=1.0_r8/rhoref

      DO iv1=1,NCS
       DO iv2=1,NCS
        DO iv3=iv2,NCS
           ! fragmentation after collision probability based on 
	   ! Gval for particles iv2 and iv3
           ! gcolfrag=(collision induced shear) / (floc strength)

           gcolfragmin=2.0_r8*(Gval*(fdiam(iv2)+fdiam(iv3)))            &
     &          **2.0_r8*fmass(iv2)*fmass(iv3)/(pi*f_fy*f_fp*           &
     &          fdiam(iv3)**2.0_r8*(fmass(iv2)+fmass(iv3))              &
     &          *((SEDFLOCS(ng)%f_rho(iv3)-rhoref)*cff2)**cff1)

           gcolfragmax=2.0_r8*(Gval*(fdiam(iv2)+fdiam(iv3)))            &
     &          **2.0_r8*fmass(iv2)*fmass(iv3)/(pi*f_fy*f_fp*           &
     &          fdiam(iv2)**2.0_r8*(fmass(iv2)+fmass(iv3))              &
     &          *((SEDFLOCS(ng)%f_rho(iv2)-rhoref)*cff2)**cff1)


           ! first case : iv3 not eroded, iv2 eroded forming 2 particles
	   ! : iv3+f_cfcst*iv2 / iv2-f_cfcst*iv2
           IF (gcolfragmin.lt.1.0_r8 .and. gcolfragmax.ge.1.0_r8) THEN  

              IF (((fmass(iv3)+f_cfcst*fmass(iv2)).gt.fmass(iv1-1))     &
     &             .and. ((fmass(iv3)+f_cfcst*fmass(iv2)).le.           &
     &             fmass(iv1))) THEN

                 f_weight=((fmass(iv3)+f_cfcst*fmass(iv2)-              &
     &                     fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (fmass(iv3)+f_cfcst*fmass(iv2).gt.fmass(iv1)       &
     &              .and. fmass(iv3)+f_cfcst*fmass(iv2).lt.             &
     &              fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_r8
                 ELSE

                    f_weight=1.0_r8-((fmass(iv3)+f_cfcst*fmass(iv2)-    &
     &                      fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_r8
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS(ng)%f_coll_prob_sh(iv2,iv3))*(fmass(iv3)+   &
     &            f_cfcst*fmass(iv2))/fmass(iv1)

              IF (fmass(iv2)-f_cfcst*fmass(iv2).gt.fmass(iv1-1)         &
     &            .and. fmass(iv2)-f_cfcst*fmass(iv2).le.               &
     &            fmass(iv1)) THEN

                 f_weight=((fmass(iv2)-f_cfcst*fmass(iv2)-              &
     &                     fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (fmass(iv2)-f_cfcst*fmass(iv2).gt.fmass(iv1)       &
     &            .and. fmass(iv2)-f_cfcst*fmass(iv2).lt.               &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_r8
                 ELSE

                    f_weight=1.0_r8-((fmass(iv2)-f_cfcst*fmass(iv2)-    &
     &                       fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_r8
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS(ng)%f_coll_prob_sh(iv2,iv3))*(fmass(iv2)-   &
     &            f_cfcst*fmass(iv2))/fmass(iv1)


              ! second case : iv3 eroded and iv2 eroded forming 3 particles : 
	      !iv3-f_cfcst*iv3 / iv2-f_cfcst*iv2 / f_cfcst*iv3+f_cfcst*iv2
           ELSEIF (gcolfragmin.ge.1.0_r8 .and. gcolfragmax.ge.          &
     &            1.0_r8) THEN  

              IF (f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3).gt.             &
     &            fmass(iv1-1) .and. f_cfcst*fmass(iv2)+f_cfcst*        &
     &            fmass(iv3).le.fmass(iv1)) THEN

                 f_weight=((f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3)-      &
     &                   fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1)))

              ELSEIF (f_cfcst*fmass(iv2)+f_cfcst*fmass(iv3).gt.         &
     &            fmass(iv1) .and. f_cfcst*fmass(iv2)+f_cfcst*          & 
     &            fmass(iv3).lt.fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_r8
                 ELSE

                    f_weight=1.0_r8-((f_cfcst*fmass(iv2)+f_cfcst*       &
     &                 fmass(iv3)-fmass(iv1))/(fmass(iv1+1)-            &
     &                 fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_r8
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &            (SEDFLOCS(ng)%f_coll_prob_sh(iv2,iv3))*(f_cfcst*      &
     &            fmass(iv2)+f_cfcst*fmass(iv3))/fmass(iv1)

              IF ((1.0_r8-f_cfcst)*fmass(iv2).gt.fmass(iv1-1)           &
     &            .and. (1.0_r8-f_cfcst)*fmass(iv2).le.                 &
     &            fmass(iv1)) THEN

                 f_weight=((1.0_r8-f_cfcst)*fmass(iv2)-                 &
     &                    fmass(iv1-1))/(fmass(iv1)-fmass(iv1-1))

              ELSEIF ((1.0_r8-f_cfcst)*fmass(iv2).gt.fmass(iv1)         &
     &            .and. (1.0_r8-f_cfcst)*fmass(iv2).lt.                 &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_r8
                 ELSE

                    f_weight=1.0_r8-(((1.0_r8-f_cfcst)*fmass(iv2)-      &
     &                      fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_r8
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &                  (SEDFLOCS(ng)%f_coll_prob_sh(iv2,iv3))*         &
     &                  ((1.0_r8-f_cfcst)*fmass(iv2))/fmass(iv1) 


              IF ((1.0_r8-f_cfcst)*fmass(iv3).gt.fmass(iv1-1) .and.     &
     &            (1.0_r8-f_cfcst)*fmass(iv3).le.fmass(iv1)) THEN

                 f_weight=((1.0_r8-f_cfcst)*fmass(iv3)-fmass(iv1-1))    &
     &                   /(fmass(iv1)-fmass(iv1-1))

              ELSEIF ((1.0_r8-f_cfcst)*fmass(iv3).gt.fmass(iv1)         &
     &            .and. (1.0_r8-f_cfcst)*fmass(iv3).lt.                 &
     &            fmass(iv1+1)) THEN

                 IF (iv1.eq.NCS) THEN 
                    f_weight=1.0_r8
                 ELSE

                    f_weight=1.0_r8-(((1.0_r8-f_cfcst)*fmass(iv3)-      &
     &                       fmass(iv1))/(fmass(iv1+1)-fmass(iv1)))
                 ENDIF

              ELSE
                 f_weight=0.0_r8
              ENDIF

              f_g4(iv2,iv3,iv1)=f_g4(iv2,iv3,iv1)+f_weight*             &
     &                  (SEDFLOCS(ng)%f_coll_prob_sh(iv2,iv3))*         &
     &                  ((1.0_r8-f_cfcst)*fmass(iv3))/fmass(iv1)


           ENDIF ! end collision test case
        ENDDO
       ENDDO
      ENDDO

      DO iv1=1,NCS
       DO iv2=1,NCS

        gcolfragiv1=2.0_r8*(Gval*(fdiam(iv1)+fdiam(iv2)))**2.0_r8*      &
     &            fmass(iv1)*fmass(iv2)/(pi*f_fy*f_fp*fdiam(iv1)        &
     &            **2.0_r8*(fmass(iv1)+fmass(iv2))*                     &
     &            ((SEDFLOCS(ng)%f_rho(iv1)-                            &
     &            rhoref)*cff2)**cff1)

        gcolfragiv2=2.0_r8*(Gval*(fdiam(iv1)+fdiam(iv2)))**2.0_r8*      &
     &            fmass(iv1)*fmass(iv2)/(pi*f_fy*f_fp*fdiam(iv2)        &
     &            **2.0_r8*(fmass(iv1)+fmass(iv2))*                     &
     &            ((SEDFLOCS(ng)%f_rho(iv2)-                            &
     &            rhoref)*cff2)**cff1)    

        mult=1.0_r8
        IF (iv1.eq.iv2) mult=2.0_r8

        IF (iv1.eq.MAX(iv1,iv2) .and. gcolfragiv1.ge.1.0_r8) THEN
           f_l4(iv2,iv1)=f_l4(iv2,iv1)+mult*                            &
     &            (SEDFLOCS(ng)%f_coll_prob_sh(iv1,iv2))
        ELSEIF (iv1.eq.MIN(iv1,iv2) .and. gcolfragiv2.ge.1.0_r8) THEN
           f_l4(iv2,iv1)=f_l4(iv2,iv1)+mult*                            &
     &            (SEDFLOCS(ng)%f_coll_prob_sh(iv1,iv2))
        ENDIF

       ENDDO
      ENDDO

      f_g4(1:NCS,1:NCS,1:NCS)=                                          &
     &            f_g4(1:NCS,1:NCS,1:NCS)*f_collfragparam
      f_l4(1:NCS,1:NCS)=f_l4(1:NCS,1:NCS)*f_collfragparam

      RETURN
      END SUBROUTINE flocmod_collfrag

!!===========================================================================
      SUBROUTINE flocmod_comp_fsd(NNin,NNout,Gval,f_dt,ng)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_comp_fsd  ***
  !&E
  !&E ** Purpose : computation of floc size distribution
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
      USE mod_param
      USE mod_scalars
      USE mod_sedflocs
!
      implicit none 

  !! * Arguments
      integer, intent(in) :: ng
      real(r8),intent(in) :: Gval
      real(r8),dimension(1:NCS),intent(in)  :: NNin
      real(r8),dimension(1:NCS),intent(out) :: NNout
      real(r8),intent(in) :: f_dt

  !! * Local declarations
      integer      :: iv1,iv2,iv3
      real(r8) :: tmp_g1,tmp_g3,tmp_l1,tmp_l3,tmp_l4,tmp_g4
      real(r8),dimension(1:NCS,1:NCS,1:NCS)     :: f_g1_tmp,f_g4
      real(r8),dimension(1:NCS,1:NCS)           :: f_l1_tmp,f_l4
      
  !!--------------------------------------------------------------------------
  !! * Executable part

      tmp_g1=0.0_r8
      tmp_g3=0.0_r8
      tmp_g4=0.0_r8
      tmp_l1=0.0_r8
      tmp_l3=0.0_r8
      tmp_l4=0.0_r8    
      f_g1_tmp(1:NCS,1:NCS,1:NCS)=0.0_r8
      f_l1_tmp(1:NCS,1:NCS)=0.0_r8

      IF (l_COLLFRAG) CALL flocmod_collfrag(Gval,ng,f_g4,f_l4)

      DO iv1=1,NCS
       DO iv2=1,NCS
        DO iv3=1,NCS
           IF (l_ASH) THEN
              f_g1_tmp(iv2,iv3,iv1)=f_g1_tmp(iv2,iv3,iv1)+              &
     &            SEDFLOCS(ng)%f_g1_sh(iv2,iv3,iv1)*Gval   
           ENDIF
           IF (l_ADS) THEN
              f_g1_tmp(iv2,iv3,iv1)=f_g1_tmp(iv2,iv3,iv1)+              &
     &            SEDFLOCS(ng)%f_g1_ds(iv2,iv3,iv1)   
           ENDIF

           tmp_g1=tmp_g1+(NNin(iv3)*(f_g1_tmp(iv2,iv3,iv1))*NNin(iv2))

           IF (l_COLLFRAG) THEN
              tmp_g4=tmp_g4+(NNin(iv3)*                                 &
     &            (SEDFLOCS(ng)%f_g4(iv2,iv3,iv1)*Gval)*NNin(iv2))
           ENDIF
        ENDDO

        tmp_g3=tmp_g3+SEDFLOCS(ng)%f_g3(iv2,iv1)*NNin(iv2)*             &
     &            Gval**1.5_r8

        IF (l_ASH) THEN
           f_l1_tmp(iv2,iv1)=f_l1_tmp(iv2,iv1)+                         &
     &            SEDFLOCS(ng)%f_l1_sh(iv2,iv1)*Gval
        ENDIF
        IF (l_ADS) THEN
           f_l1_tmp(iv2,iv1)=f_l1_tmp(iv2,iv1)+                         &
     &            SEDFLOCS(ng)%f_l1_ds(iv2,iv1)*Gval
        ENDIF

        tmp_l1=tmp_l1+(f_l1_tmp(iv2,iv1))*NNin(iv2)

        IF (l_COLLFRAG) THEN
           tmp_l4=tmp_l4+(SEDFLOCS(ng)%f_l4(iv2,iv1)*Gval)*NNin(iv2)
        ENDIF
       ENDDO

       tmp_l1=tmp_l1*NNin(iv1)
       tmp_l4=tmp_l4*NNin(iv1)
       tmp_l3=SEDFLOCS(ng)%f_l3(iv1)*Gval**1.5_r8*NNin(iv1)

       NNout(iv1)=NNin(iv1)+ f_dt*(tmp_g1+tmp_g3+tmp_g4-(tmp_l1+        &
     &            tmp_l3+tmp_l4))

       tmp_g1=0.0_r8
       tmp_g3=0.0_r8
       tmp_g4=0.0_r8
       tmp_l1=0.0_r8
       tmp_l3=0.0_r8
       tmp_l4=0.0_r8    
      ENDDO

      RETURN
      END SUBROUTINE flocmod_comp_fsd



!!===========================================================================
      SUBROUTINE flocmod_mass_control(NN,mneg,ng)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_mass_control  ***
  !&E
  !&E ** Purpose : Compute mass in every class after flocculation and 
  !&E              returns negative mass if any
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
      USE mod_sedflocs
      USE mod_param
      USE mod_scalars
!
      implicit none 
      integer, intent(in) :: ng

  !! * Local declarations
      integer      :: iv1
      real(r8),intent(out)     :: mneg
      real(r8),dimension(1:NCS),intent(in)     :: NN
      !real(r8),DIMENSION(0:NCS+1)      :: f_mass

  !!--------------------------------------------------------------------------
  !! * Executable part

      mneg=0.0_r8

      DO iv1=1,NCS
       IF (NN(iv1).lt.0.0_r8) THEN
         mneg=mneg-NN(iv1)*SEDFLOCS(ng)%f_mass(iv1)
       ENDIF
      ENDDO

      RETURN
      END SUBROUTINE flocmod_mass_control

!!===========================================================================
      SUBROUTINE flocmod_mass_redistribute(NN,ng)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_mass_redistribute  ***
  !&E
  !&E ** Purpose : based on a tolerated negative mass parameter, negative masses  
  !&E              are redistributed linearly towards remaining postive masses 
  !&E              and negative masses are set to 0
  !&E                   
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
      USE mod_param
      USE mod_scalars
      USE mod_sedflocs
!
      implicit none 
      integer, intent(in) :: ng


  !! * Local declarations
      integer      :: iv
      real(r8)     :: npos
      real(r8)     :: mneg
      real(r8),dimension(1:NCS),intent(inout)     :: NN
      real(r8),dimension(1:NCS)                   :: NNtmp
      !real(r8),DIMENSION(0:NCS+1)      :: f_mass
  !!--------------------------------------------------------------------------
  !! * Executable part

      mneg=0.0_r8
      npos=0.0_r8
      NNtmp(:)=NN(:)

      DO iv=1,NCS
       IF (NN(iv).lt.0.0_r8) THEN
         mneg=mneg-NN(iv)*SEDFLOCS(ng)%f_mass(iv)
         NNtmp(iv)=0.0_r8
       ELSE
         npos=npos+1.0_r8
       ENDIF
      ENDDO

      IF (mneg.gt.0.0_r8) THEN 
       IF (npos.eq.0.0_r8) THEN
         WRITE(*,*) 'CAUTION : all floc sizes have negative mass!'
         WRITE(*,*) 'SIMULATION STOPPED'
         STOP    
       ELSE
         DO iv=1,NCS
           IF (NN(iv).gt.0.0_r8) THEN
              NN(iv)=NN(iv)-mneg/sum(NNtmp)*NN(iv)/                     &
     &            SEDFLOCS(ng)%f_mass(iv)
           ELSE
              NN(iv)=0.0_r8
           ENDIF

         ENDDO

       ENDIF
      ENDIF

      RETURN
      END SUBROUTINE flocmod_mass_redistribute

!!===========================================================================    
      SUBROUTINE flocmod_comp_g(k,i,j,Gval,diss,ng)

  !&E--------------------------------------------------------------------------
  !&E                 ***  ROUTINE flocmod_comp_g  ***
  !&E
  !&E ** Purpose : compute shear rate to estimate shear aggregation and erosion  
  !&E 
  !&E ** Description :
  !&E
  !&E ** Called by : flocmod_main
  !&E
  !&E ** External calls : 
  !&E
  !&E ** Reference :
  !&E
  !&E ** History :
  !&E     ! 2013-09 (Romaric Verney)
  !&E
  !&E--------------------------------------------------------------------------
  !! * Modules used
      USE mod_sedflocs
      USE mod_param
      USE mod_scalars
!
      implicit none 

  !! * Local declarations
      integer,  intent(in)      :: k,i,j
      integer, intent(in) :: ng
      real(r8),intent(out)     :: Gval
      real(r8)     :: htn,ustar,z,diss,nueau
! l_testcase - if .TRUE. sets G(t) to values from lab experiment
!      logical, parameter :: l_testcase = .TRUE.
  !!--------------------------------------------------------------------------
  !! * Executable part
  !    nueau=1.06e-6_r8
  !    ustar=sqrt(tenfon(i,j)/rhoref)
  !    htn=h0(i,j)+ssh(i,j)
  !    z=(1.0_r8+sig(k))*htn
      
!
! ALA from CRS
!
       nueau = 1.5E-6_r8
      IF (l_testcase) THEN
        ! reproducing flocculation experiment Verney et al., 2011
       Gval=0.0_r8
       IF (time(ng) .lt. 7201.0_r8) THEN
        Gval=1.0_r8
       ELSEIF (time(ng) .lt. 8401.0_r8) THEN
        Gval=2.0_r8
       ELSEIF (time(ng) .lt. 9601.0_r8) THEN
        Gval=3.0_r8  
       ELSEIF (time(ng) .lt. 10801.0_r8) THEN
        Gval=4.0_r8
       ELSEIF (time(ng) .lt. 12601.0_r8) THEN
        Gval=12.0_r8
       ELSEIF (time(ng) .lt. 13801.0_r8) THEN
        Gval=4.0_r8
       ELSEIF (time(ng) .lt. 15001.0_r8) THEN
        Gval=3.0_r8
       ELSEIF (time(ng) .lt. 16201.0_r8) THEN
        Gval=2.0_r8
       ELSEIF (time(ng) .lt. 21601.0_r8) THEN
        Gval=1.0_r8
       ELSEIF (time(ng) .lt. 25201.0_r8) THEN
        Gval=0.0_r8
       ELSEIF (time(ng) .lt. 30601.0_r8) THEN
        Gval=1.0_r8
       ELSEIF (time(ng) .lt. 31801.0_r8) THEN
        Gval=2.0_r8                     
       ELSEIF (time(ng) .lt. 33001.0_r8) THEN
        Gval=3.0_r8       
       ELSEIF (time(ng) .lt. 34201.0_r8) THEN
        Gval=4.0_r8
       ELSEIF (time(ng) .lt. 36001.0_r8) THEN
        Gval=12.0_r8
       ELSEIF (time(ng) .lt. 37201.0_r8) THEN
        Gval=4.0_r8
       ELSEIF (time(ng) .lt. 38401.0_r8) THEN
        Gval=3.0_r8
       ELSEIF (time(ng) .lt. 39601.0_r8) THEN
        Gval=2.0_r8
       ELSEIF (time(ng) .lt. 45001.0_r8) THEN
        Gval=1.0_r8
       ELSEIF (time(ng) .lt. 48601.0_r8) THEN
        Gval=0.0_r8
       ELSEIF (time(ng) .lt. 54001.0_r8) THEN
        Gval=1.0_r8
       ELSEIF (time(ng) .lt. 55201.0_r8) THEN
        Gval=2.0_r8                     
       ELSEIF (time(ng) .lt. 56401.0_r8) THEN
        Gval=3.0_r8       
       ELSEIF (time(ng) .lt. 57601.0_r8) THEN
        Gval=4.0_r8
       ELSE 
        Gval=12.0_r8
       ENDIF
      ELSE
         Gval=sqrt(diss/nueau)
      ENDIF
! NO KLUDGE
!      Gval = 12.0_r8
      RETURN
      END SUBROUTINE flocmod_comp_g
      
      
#endif
      END MODULE sed_flocs_mod