init_transport.f

Go to the documentation of this file.

!                   *************************
                    SUBROUTINE init_transport
!                   *************************
!
     &(trouve,debu,hiding,nsicla,npoin,
     & t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,
     & charr,qs_c,qsxc,qsyc,calfa_cl,salfa_cl,coefpn,slopeff,
     & susp,qs_s,qs,qscl,qscl_c,qscl_s,qsclxs,qsclys,
     & unorm,u2d,v2d,hn,cf,mu,tob,tobw,uw,tw,thetaw,fw,houle,
     & avail,acladm,unladm,ksp,ksr,ks,
     & icf,hidfac,xmvs,xmve,grav,vce,hmin,karman,
     & zero,pi,ac,imp_inflow_c,zref,icq,cstaeq,csratio,
     & cmax,cs,cs0,seccurrent,bijk,
     & ielmt,fdm,xwc,fd90,sedco,vitce,partheniades,vitcd,
     & u3d,v3d,code)
!
!***********************************************************************
! SISYPHE   V7P2                                   21/07/2011
!***********************************************************************
!
!brief
!
!history  JMH
!+        24/01/2008
!+
!+   TEST 'IF(CHARR.OR.SUSP)' ADDED AT THE END
!
!history  JMH
!+        16/09/2009
!+
!+   AVAIL(NPOIN,10,NSICLA)
!
!history  JMH
!+        07/12/2009
!+
!+   MODIFICATIONS FOR RESTART WITH WARNINGS WHEN A VARIABLE
!
!history  C. VILLARET (LNHE)
!+
!+        V6P0
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!+
!history  C. VILLARET (LNHE)+ R. KOPMANN + U.MERKEL
!+
!+        20/03/2011
!+        V6P1
!
!history  C.VILLARET (EDF-LNHE), P.TASSI (EDF-LNHE)
!+        19/07/2011
!+        V6P1
!+  Name of variables
!
!history  MAK (HRW)
!+        31/05/2012
!+        V6P2
!+  Include CSRATIO
!
!history  P TASSI (LNHE)
!+        18/06/2012
!+        V6P2
!+   updated version with HRW's development for Soulsby-van Rijn's concentration
!
!history  R KOPMANN (BAW)
!+        10/05/2016
!+        V7P2
!+ CALFA,SALFA dependent of grain classes
!
!history  F.CORDIER & P.TASSI (EDF-LNHE)
!+        12/09/2018
!+        V8P0
!+  CALCULATION OF SAND FRACTION CONTENT AT EACH NODE
!+  FOR WILCOCK AND CROWE FORMULA (2003)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| CRITICAL SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| AVAIL          |<->| VOLUME PERCENT OF EACH CLASS
!| BIJK           |-->| COEFFICIENT OF THE BIJKER FORMULA
!| CALFA          |<->| COSINUS OF THE ANGLE BETWEEN MEAN FLOW AND TRANSPORT
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| CHARR          |-->| BEDLOAD
!| CMAX           |---| MAX(PARTHENIADES/SETTLING VELOCITY)
!| COEFPN         |<->| CORRECTION OF TRANSORT FOR SLOPING BED EFFECT
!| CS             |<->| CONCENTRATION AT TIME N
!| CS0            |-->| CONCENTRATION AT TIME 0
!| CSTAEQ         |<->| EQUILIBRIUM CONCENTRATION
!| CSRATIO        |<->| EQUILIBRIUM CONCENTRATION FOR SOULSBY-VAN RIJN EQ.
!| DEBUG          |-->| FLAG FOR DEBUGGING
!| FD90           |-->| DIAMETER D90
!| FDM            |-->| DIAMETER DM FOR EACH CLASS
!| FW             |-->| WAVE FRICTION FACTOR
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HIDFAC         |-->| HIDING FACTOR FORMULAS
!| HIDING         |-->| HIDING FACTOR CORRECTION
!| HMIN           |-->| MINIMUM VALUE OF WATER DEPTH
!| HN             |-->| WATER DEPTH
!| HOULE          |-->| LOGICAL, FOR WAVE EFFECTS
!| ICF            |-->| BED-LOAD OR TOTAL LOAD TRANSPORT FORMULAS
!| ICQ            |-->| REFERENCE CONCENTRATION FORMULA
!| IELMT          |-->| NUMBER OF ELEMENTS
!| IMP_INFLOW_C   |-->| IMPOSED CONCENTRATION IN INFLOW
!| KARMAN         |-->| VON KARMAN CONSTANT
!| KS             |-->| BED ROUGHNESS
!| KSP            |-->| BED SKIN ROUGHNESS
!| KSR            |-->| RIPPLE BED ROUGHNESS
!| MU             |<->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NPOIN          |-->| NUMBER OF POINTS
!| NSICLA         |-->| NUMBER OF SEDIMENT CLASSES
!| PARTHENIADES   |-->| CONSTANT OF THE KRONE AND PARTHENIADES EROSION LAW (KG/M2/S)
!| PI             |-->| PI
!| QS             |<->| BEDLOAD TRANSPORT RATE
!| QSCL           |<->| SUSPENDED LOAD TRANSPORT RATE
!| QSCLXS         |<->| SUSPENDED LOAD TRANSPORT RATE FOR EACH CLASS X-DIRECTION
!| QSCLYS         |<->| SUSPENDED LOAD TRANSPORT RATE FOR EACH CLASS Y-DIRECTION
!| QSCL_C         |<->| BEDLOAD TRANSPORT RATE
!| QSCL_S         |<->| SUSPENDED LOAD TRANSPORT RATE
!| QSXS           |<->| SOLID DISCHARGE X (SUSPENSION)
!| QSYS           |<->| SOLID DISCHARGE Y (SUSPENSION)
!| QS_C           |-->| BEDLOAD TRANSPORT RATE
!| QS_S           |<->| SUSPENDED LOAD TRANSPORT RATE
!| SALFA          |<->| SINUS OF THE ANGLE BETWEEN TRANSPORT RATE AND CURRENT
!| SECCURRENT     |-->| LOGICAL, PARAMETRISATION FOR SECONDARY CURRENTS
!| SEDCO          |-->| LOGICAL, SEDIMENT COHESIVE OR NOT
!| SLOPEFF        |-->| LOGICAL, SLOPING BED EFFECT OR NOT
!| SUSP           |-->| LOGICAL, SUSPENSION
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T10            |<->| WORK BIEF_OBJ STRUCTURE
!| T11            |<->| WORK BIEF_OBJ STRUCTURE
!| T12            |<->| WORK BIEF_OBJ STRUCTURE
!| T13            |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| T4             |<->| WORK BIEF_OBJ STRUCTURE
!| T5             |<->| WORK BIEF_OBJ STRUCTURE
!| T6             |<->| WORK BIEF_OBJ STRUCTURE
!| T7             |<->| WORK BIEF_OBJ STRUCTURE
!| T8             |<->| WORK BIEF_OBJ STRUCTURE
!| T9             |<->| WORK BIEF_OBJ STRUCTURE
!| THETAW         |-->| ANGLE BETWEEN WAVE AND CURRENT
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOBW           |-->| WAVE INDUCED SHEAR STRESS
!| TW             |-->| WAVE PERIOD
!| U2D            |<->| MEAN FLOW VELOCITY X-DIRECTION
!| UNLADM         |-->| MEAN DIAMETER OF ACTIVE STRATUM LAYER
!| UNORM          |<->| NORM OF THE MEAN FLOW VELOCITY
!| UW             |-->| ORBITAL WAVE VELOCITY
!| V2D            |<->| MEAN FLOW VELOCITY Y-DIRECTION
!| VCE            |-->| WATER VISCOSITY
!| VITCD          |-->| CRITICAL SHEAR VELOCITY FOR MUD DEPOSITION
!| VITCE          |-->| CRITICAL EROSION SHEAR VELOCITY OF THE MUD
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| WATER DENSITY
!| XWC            |-->| SETTLING VELOCITY
!| ZERO           |-->| ZERO
!| ZREF           |-->| REFERENCE ELEVATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE interface_sisyphe, ex_init_transport => init_transport
!
      USE declarations_sisyphe, ONLY : nomblay,mpm_aray,mpm
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)              :: NSICLA,NPOIN,TROUVE(*),ICQ
      INTEGER, INTENT(IN)              :: ICF,HIDFAC,IELMT,SLOPEFF
      LOGICAL, INTENT(IN)              :: CHARR,DEBU,SUSP,IMP_INFLOW_C
      LOGICAL, INTENT(IN)              :: SECCURRENT,SEDCO(*)
      LOGICAL, INTENT(IN)              :: HOULE
      TYPE(bief_obj),    INTENT(IN)    :: U2D,V2D,UNORM,HN,CF
      TYPE(bief_obj),    INTENT(IN)    :: MU,TOB,TOBW,UW,TW,THETAW,FW
      TYPE(bief_obj),    INTENT(IN)    :: ACLADM,UNLADM,KSP,KSR,KS
      TYPE(bief_obj),    INTENT(INOUT) :: HIDING
      TYPE(bief_obj),    INTENT(INOUT) :: QS_C, QSXC, QSYC
      TYPE(bief_obj),    INTENT(INOUT) ::  CALFA_CL,SALFA_CL
      TYPE(bief_obj),    INTENT(INOUT) :: T1,T2,T3,T4,T5,T6,T7,T8
      TYPE(bief_obj),    INTENT(INOUT) :: T9,T10,T11,T12
      TYPE(bief_obj),    INTENT(INOUT) :: ZREF,CSTAEQ,CSRATIO
      TYPE(bief_obj),    INTENT(INOUT) :: CS
      TYPE(bief_obj),    INTENT(INOUT) :: QS_S,QS,QSCL_C,QSCL_S
      TYPE(bief_obj),    INTENT(INOUT) :: COEFPN
      TYPE(bief_obj),    INTENT(INOUT) :: QSCLXS,QSCLYS,QSCL
      DOUBLE PRECISION,  INTENT(IN)    :: XMVS,XMVE,GRAV,VCE
      DOUBLE PRECISION,  INTENT(IN)    :: HMIN,KARMAN,ZERO,PI
      DOUBLE PRECISION,  INTENT(IN)    :: PARTHENIADES,BIJK,XWC(nsicla)
      DOUBLE PRECISION,  INTENT(IN)    :: FD90(nsicla),CS0(nsicla)
      DOUBLE PRECISION,  INTENT(IN)    :: VITCE,VITCD
      DOUBLE PRECISION,  INTENT(INOUT) :: AC(nsicla),CMAX,FDM(nsicla)
      DOUBLE PRECISION,  INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
!
      TYPE(bief_obj),    INTENT(IN)    :: U3D,V3D
      CHARACTER(LEN=24), INTENT(IN)    :: CODE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,J
      DOUBLE PRECISION AAA,USTARP,U3DNORM
      LOGICAL NEED_CS
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION SANFRA(npoin)
      INTEGER K
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
! --- START : INITIALISES RATE OF TRANSPORT AND SUSPENSION
!
!     FOR INITIALISATION : SLOPE EFFECT AND DEVIATION ARE CANCELLED
!
!     RK in case of coupling with T3D, the direction should
!     come from the bottom velocity
!
!     Calculation of sand fraction content at each node (Wilcock and Crowe, 2003)
      IF (icf == 10) THEN
        DO k = 1, npoin
          sanfra(k) = 0.0d0
          DO i = 1, nsicla
            IF (fdm(i).LT.2d-3) THEN
              sanfra(k) = sanfra(k) + avail(k,1,i)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
!
      IF(code(1:9).EQ.'TELEMAC3D') THEN
        DO i=1,npoin
          u3dnorm=sqrt(u3d%R(i)*u3d%R(i)+v3d%R(i)*v3d%R(i))
          IF(u3dnorm.GE.1.d-12) THEN
            calfa_cl%ADR(1)%P%R(i)=u3d%R(i)/u3dnorm
            salfa_cl%ADR(1)%P%R(i)=v3d%R(i)/u3dnorm
          ELSE
            calfa_cl%ADR(1)%P%R(i)=1.d0
            salfa_cl%ADR(1)%P%R(i)=0.d0
          ENDIF
        ENDDO
      ELSE
        CALL os('X=Y/Z   ',x=calfa_cl%ADR(1)%P, y=u2d, z=unorm,
     &           c=0.d0, iopt=2, infini=1.d0, zero=1.d-12)
        CALL os('X=Y/Z   ',x=salfa_cl%ADR(1)%P, y=v2d, z=unorm,
     &           c=0.d0, iopt=2, infini=0.d0, zero=1.d-12)
      ENDIF
      IF(nsicla.GT.1) THEN
        DO i=2,nsicla
          CALL os('X=Y     ', x=calfa_cl%ADR(i)%P,
     &            y=calfa_cl%ADR(1)%P)
          CALL os('X=Y     ', x=salfa_cl%ADR(i)%P,
     &            y=salfa_cl%ADR(1)%P)
        ENDDO
      ENDIF
!
!     appel a effpnt ?
!
      CALL os('X=C     ',x=coefpn,c=1.d0)
!
      IF(charr) THEN
!
!       MPM for each Layer
!
        CALL os('X=C     ', x=mpm_aray, c=mpm)
!
        CALL os('X=C     ',x=hiding,c=1.d0)
!
        DO i = 1, nsicla
!
          IF(sedco(i)) THEN
!           IF COHESIVE: NO BEDLOAD TRANSPORT
            CALL os('X=0     ', x=qscl_c%ADR(i)%P)
          ELSE
!           IF NON COHESIVE
            CALL bedload_formula
     &        (u2d,v2d,unorm,hn,cf,mu,tob,tobw,uw,tw,thetaw,fw,
     &        acladm, unladm,ksp,ksr,avail(1:npoin,1,i),
     &        npoin,icf,hidfac,xmvs,xmve,
     &        fdm(i),grav,vce,hmin,xwc(i),fd90(i),karman,zero,
     &        pi,susp,ac(i),hiding,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
     &        t11,t12,qscl_c%ADR(i)%P,qscl_s%ADR(i)%P,
     &        ielmt,seccurrent,slopeff,
     &        coefpn,calfa_cl%ADR(i)%P,salfa_cl%ADR(i)%P,
! FC
     &        bijk,houle,sanfra)
! ENDFC
!
          ENDIF
!         SUM ON ALL CLASSES
          DO j=1,npoin
            qs_c%R(j) = qs_c%R(j) + qscl_c%ADR(i)%P%R(j)
!
!       COMPUTES THE X AND Y COMPONENTS OF TRANSPORT
            qsxc%R(j) = qsxc%R(j) + qscl_c%ADR(i)%P%R(j)
     &                  *calfa_cl%ADR(i)%P%R(j)
            qsyc%R(j) = qsyc%R(j) + qscl_c%ADR(i)%P%R(j)
     &                  *salfa_cl%ADR(i)%P%R(j)
          ENDDO
!
        ENDDO
!
!
      ENDIF
!
!     START : COMPUTES THE SUSPENDED TRANSPORT
!
      IF(susp) THEN
!
!       INITIALISES ZREF
!
        IF(icq.EQ.1) THEN
          CALL os('X=Y     ', x=zref, y=ksp)
        ELSEIF(icq.EQ.2) THEN
          CALL os('X=Y     ', x=zref, y=ksr)
        ELSEIF(icq.EQ.3) THEN
          CALL os('X=CY    ', x=zref, y=ks,c=0.5d0)
        ENDIF
!
!       FOR RANK OF CS IN TROUVE SEE POINT_SISYPHE, NOMVAR_SISYPHE
!       22+I+(NOMBLAY+1)*NSICLA IS THE ADDRESS OF CONCENTRATIONS
!
        need_cs=.false.
        DO i=1,nsicla
          IF(trouve(22+i+(nomblay+1)*nsicla).EQ.0) need_cs=.true.
        ENDDO
!
!       COMPUTES THE INITIAL CONCENTRATIONS
!
        IF(.NOT.debu.OR.need_cs) THEN
!
          CALL condim_susp(cs,cs0,nsicla)
!
!         END MODIFICATIONS (CV)
!         OPTION: IMPOSED INFLOW CONCENTRATIONS ...
!
          IF(imp_inflow_c) THEN
!
!           TAUP IN T8
            CALL os('X=CYZ   ', x=t8, y=tob, z=mu, c=1.d0)
!           USTAR (TOTAL) IN T9
            CALL os('X=CY    ', x=t9, y=tob, c=1.d0/xmve)
            CALL os('X=SQR(Y)', x=t9, y=t9)
!
!           START: LOOP ON THE CLASSES
!
            DO i=1,nsicla
!
              IF(.NOT.sedco(i)) THEN
!
!               NON COHESIVE SED: INITIALISES CSTAEQ
!
                IF(icq.EQ.1) THEN
                  CALL suspension_fredsoe( fdm(i) ,t8,npoin,
     &                grav, xmve, xmvs, ac(i),  cstaeq )
                ELSEIF(icq.EQ.2) THEN
                  CALL suspension_bijker(t8,npoin,charr,qs_c,zref,
     &                                   zero,cstaeq,xmve)
                ELSEIF(icq.EQ.3) THEN
                  CALL suspension_vanrijn(fdm(i),t8,npoin,
     &               grav,xmve,xmvs,vce,zero,ac(i), cstaeq,zref)
                ELSEIF(icq.EQ.4) THEN
                  csratio%R=1d0
                  CALL suspension_sandflow(fdm(i),fd90(i),npoin,
     &                                     grav,xmve,xmvs,zero,
     &                                     cstaeq,hn,u2d,v2d,
     &                                     csratio)
                ENDIF
!               ROUSE CONCENTRATION PROFILE IS ASSUMED BASED ON TOTAL FRICTION
!               VELOCITY
!
                CALL suspension_rouse(t9,hn,npoin,
     &                             karman,zero,xwc(i),zref,t12)
!
                DO j=1,npoin
                  cstaeq%R(j)=cstaeq%R(j)*avail(j,1,i)
                ENDDO
!               CALL OS( 'X=XY    ',X=CSTAEQ,Y=AVAI%ADR(I)%P)
                CALL os( 'X=Y/Z   ',x=cs%ADR(i)%P,y=cstaeq,z=t12)
!
!               END NON-COHESIVE
!
              ELSE
!
!               FOR COHESIVE SEDIMENT
!
!               THIS VALUE CAN BE ALSO CHANGED BY THE USER
!               IN SUBROUTINE USER_KRONE_PARTHENIADES
!
                CALL os('X=Y     ', x=zref, y=ksp)
!
                cmax = max(cmax,partheniades/xwc(i))
!
                IF(vitce.GT.1.d-8.AND.vitcd.GT.1.d-8) THEN
                  DO j = 1, npoin
!                 FLUER
                  ustarp= sqrt(t8%R(j)/xmve)
                  aaa= partheniades*
     &                max(((ustarp/vitce)**2-1.d0),zero)
!                 FLUDPT
!                 BBB=XWC(I)*MAX((1.D0-(USTARP/VITCD)**2),ZERO)
!                 IF NO DEPOSITION, THE EQUILIBRIUM CONCENTRATION IS INFINITE!
                  cs%ADR(i)%P%R(j) = aaa/xwc(i)
!
                  ENDDO
                ELSE
                  CALL os('X=0     ',x=cs%ADR(i)%P)
                ENDIF
!
                DO j=1,npoin
                  cs%ADR(i)%P%R(j)=cs%ADR(i)%P%R(j)*avail(j,1,i)
                ENDDO
!
! END COHESIVE
!
              ENDIF
!
! END OF LOOP ON THE CLASSES
!
            ENDDO
!
! END OF OPTION: IMPOSED INFLOW CONCENTRATION
!
          ENDIF
!
! END OF IF(.NOT.DEBU.OR.NEED_CS.EQ.0)) THEN
!
        ENDIF
!
! COMPUTES SUSPENDED TRANSPORT
!
        DO i=1,nsicla
!                                    UCONV DONE IN SUSPENSION_COMPUTATION
!                                    HERE WE USE
!                                    U2D AS TENTATIVE VALUE OF UCONV
          CALL os('X=YZ    ',x=t11,y=u2d, z=hn)
          CALL os('X=YZ    ',x=t12,y=v2d, z=hn)
!
          CALL os('X=YZ    ',x=qsclxs%ADR(i)%P,y=cs%ADR(i)%P,z=t11)
          CALL os('X=YZ    ',x=qsclys%ADR(i)%P,y=cs%ADR(i)%P,z=t12)
!
          CALL os('X=N(Y,Z) ',x=qscl_s%ADR(i)%P,
     &                        y=qsclxs%ADR(i)%P,z=qsclys%ADR(i)%P)
!
        ENDDO

          DO j=1,npoin
            DO i=1,nsicla
            qs_s%R(j) = qs_s%R(j) + qscl_s%ADR(i)%P%R(j)
            ENDDO
          ENDDO
      ENDIF
!
! END OF SUSPENSION
!
!
!     COMPUTES THE TRANSPORT FOR EACH CLASS (IF NOT RESTART OR IF
!                                              DATA NOT FOUND)
      DO i=1, nsicla
        WRITE(lu,*) 'QSCL REINITIALISED IN INIT_TRANSPORT'
        WRITE(lu,*) 'FOR CLASS ',i
        IF(charr.AND.susp) THEN
          CALL os('X=Y+Z   ', x=qscl%ADR(i)%P,
     &            y=qscl_s%ADR(i)%P, z=qscl_c%ADR(i)%P)
        ELSEIF(charr) THEN
          CALL os('X=Y     ',x=qscl%ADR(i)%P,y=qscl_c%ADR(i)%P)
        ELSEIF(susp) THEN
          CALL os('X=Y     ',x=qscl%ADR(i)%P,y=qscl_s%ADR(i)%P)
        ENDIF
      ENDDO
!
!     COMPUTES TOTAL TRANSPORT QS
!
      WRITE(lu,*) 'QS REINITIALISED IN INIT_TRANSPORT'
      IF(charr.AND.susp) THEN
        CALL os('X=Y+Z   ',x=qs,y=qs_c,z=qs_s)
      ELSEIF(charr) THEN
        CALL os('X=Y     ',x=qs,y=qs_c)
      ELSEIF(susp) THEN
        CALL os('X=Y     ',x=qs,y=qs_s)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END