suspension_miles.f

Go to the documentation of this file.

!                   ***************************
                    SUBROUTINE suspension_miles
!                   ***************************
!
     &(hn,npoin,hmin,fdm,fd90,xwc,csratio)
!
!***********************************************************************
! SISYPHE   V6P3                                   16/05/2011
!***********************************************************************
!
!brief    COMPUTES THE BED EXCHANGE FACTOR BETA BASED ON MILES (1986)
!         FOUND IN HR WALLINGFORD REPORT: SR75.
!
!history  D. M. KELLY (HRW)
!+        29/06/2011
!+        V6P3
!+    Fisrt version added to Sisyphe library
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| HMIN           |-->|
!| HN             |-->|
!| NPOIN          |-->|
!| CSRATIO        |---| RATIO BETWEEN BOTTOM CONC. AND AVERAGE CONC.
!| XWC            |-->| FALL VELOCITY W_S
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_sisyphe, ONLY: sfon,u2d,v2d,dt
      USE interface_sisyphe,ex_suspension_miles => suspension_miles
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj),   INTENT(IN)    :: HN
      INTEGER,          INTENT(IN)    :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: FDM,FD90,XWC,HMIN
      TYPE(bief_obj),   INTENT(INOUT) :: CSRATIO
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          :: I
      DOUBLE PRECISION :: UB,USTARS,RB,BETAS
      DOUBLE PRECISION :: SFBETA(npoin),AUX,RRTPI
      DOUBLE PRECISION :: FVINV
      DOUBLE PRECISION :: DZ,TAU,TAU_SQ,UCR
!
      DOUBLE PRECISION BIEF_ERF
      EXTERNAL         bief_erf
!
      INTRINSIC sqrt,exp
!     INTRINSIC ERF
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!     COMPUTE CONSTANTS:
!
      DOUBLE PRECISION PI
      pi = 4.d0 * atan( 1.d0 )
      rrtpi = 1.d0/sqrt(pi)
      fvinv = 1.d0/xwc**2
!
!     DMK MOD 03/05/2011
!
!     MODIFYING SUSPENDED LOAD TO BE LIKE SANDFLOW AND ALSO TO TAKE
!     LAG INTO CONSIDERATION VIA A BETA COMPUTATION.
!
      DO i=1,npoin
!
        IF(fdm.LT.500.e-6) THEN
!         UCR = 0.19D0*(D50**0.1D0)*LOG10(4.D0*(HN%R(I)/FDM))    ! SANDFLOW
          ucr = 0.19d0*(fdm**0.1d0)*log10(4.d0*hn%R(i)/fd90)  ! CORRECT
        ELSE
!         UCR = 8.5D0*(D50**0.6D0)*LOG10(4.D0*(HN%R(I)/D50))      ! SANDFLOW
          ucr = 8.5d0*(fdm**0.6d0)*log10(4.d0*hn%R(i)/fd90)    ! CORRECT
        ENDIF
!
        ub = sqrt(u2d%R(i)**2+v2d%R(i)**2)
!
        IF(hn%R(i).GT.hmin.AND.ub.GT.ucr) THEN ! STOP PROBLEMS WITH LOW VELOCITIES
          ustars = 1.3d0*ub*sqrt(sfon/8.d0)
          rb     = xwc*15.d0/ustars
!         THIS IS DIFFERENT BETWEEN MILES’ PAPERS DUE TO TYPO IN DEF OF R
          betas  = rb/(1.d0-exp(-rb))
!         ACTUAL EQUATION FROM MILES (1981)  [1]
          dz     = (1.d0/6.d0)*0.4d0*hn%R(i)*ustars
!         USE WITH [1]
          tau    = xwc*sqrt(dt/(4.d0*dz))
          tau_sq = tau**2
!         EQ. (27) NOTE THIS IS ALREADY INTEGRATED WITH RESPECT TO TIME (DTS)
!
!         IF COMPILER ALOWS...
!         AUX=ERF(TAU)
!         AND IF NOT...
          aux=bief_erf(tau)
!
          sfbeta(i)=
     &      dz*fvinv*betas*(4.d0*tau_sq*(1.d0+tau_sq)*(1.d0-aux)
     &    + aux - 2.d0*tau*(1.d0+2.d0*tau_sq)*exp(-tau_sq)*rrtpi)
!
          IF(sfbeta(i)*xwc/hn%R(i).GT.1.d0) sfbeta(i)=1.d0/(xwc*hn%R(i))
!         DIVIDE BACK THROUGH BY DT AS WE WILL INTEGRATE UP  WRTT LATER
          csratio%R(i)=sfbeta(i)/dt ! FOR USE WITH EQ. (27)
!         NOTE WE RECORD BETA_S (PROFILE PARAMETER) I.E. THE RATIO OF REF LEVEL CONC
!         TO DEPTH AVERAGED CONC THIS WILL THEN BE USED IN SUSPENSION_SANDFLOW.F:
!         NB: STORED IN T14
        ELSE
          sfbeta(i) = 1.d0
          csratio%R(i) = 1.d0
        ENDIF
!
      ENDDO
!
!     END BETA FACTOR COMPUTATION
!
!======================================================================!
!======================================================================!
!
      RETURN
      END