calcs3d_micropol.f

Go to the documentation of this file.

!                    ***************************
                      SUBROUTINE calcs3d_micropol
!                    **************************
     & (npoin2,tn,texp,timp,zprop,cf,un,vn,
     &  t2_1,t2_2,t2_3,t3_1,t3_2,debug)
!
!
!***********************************************************************
! TELEMAC2D   V7P2                                        21/05/2016
!***********************************************************************
!
!brief    COMPUTES SOURCE TERMS FOR MICROPOL WAQ PROCESS IN 3D
!          WAQ PROCESS OF CODE_TRACER (MASCARET SYSTEM)
!
!history  R. ATA
!+        21/05/2016
!+        V7P0
!+       REAL CREATION
!
!history  S.E. BOURBAN (HRW)
!+        07/06/2017
!+        V7P3
!+        Indexing tracer (IND_*) to avoid conflicting naming convention
!+        between user defined tracers, water quality processes and
!+        ice processes. Introduction of the array RANK_*.
!
!history  S.E. BOURBAN (HRW)
!+        25/09/2017
!+        V7P3
!+        TEXP and TIMP are now additive to account for a variety of
!+        of sources / sinks on a given TRACER
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CCSEDIM        |-->| CONSTANT OF EXPONENTIAL DESINTEGRATION
!| CDISTRIB       |-->| COEFFICIENT OF DISTRIBUTION (KD)
!| DEBUG          |-->| IF NE.0 THEN DEBUG MODE
!| DT             |-->| TIME STEP
!| ERO            |-->| EROSION RATE
!| KDESORP        |-->| KINETIC CONSTANT OF  DESORPTION
!| NPOIN          |-->| TOTAL NUMBER OF MESH NODES
!| TAUB           |-->| BED SHEAR
!| TAUS           |-->| CRITICAL STRESS OF RESUSPENSION
!| TAUR           |-->| SEDIMENTATION CRITICAL STRESS
!| TEXP           |<--| EXPLICIT SOURCE TERMS OF TRACERS
!| TIMP           |<--| IMPLICIT SOURCE TERMS OF TRACERS
!| TN             |-->| TRACERS
!| VITCHU         |-->| SEDIMENT SETTLING VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_special
      USE interface_parallel
      USE declarations_waqtel,ONLY:ero,taur,taus,vitchu,cdistrib,
     &  ro0,kdesorp,ccsedim,
     &  ind_ss,ind_sf,ind_c,ind_css,ind_csf
      USE interface_waqtel, ex_calcs3d_micropol => calcs3d_micropol
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          , INTENT(IN   ) :: npoin2
      INTEGER          , INTENT(IN   ) :: debug
      TYPE(bief_obj)   , INTENT(IN   ) :: tn,zprop,cf,un,vn
      TYPE(bief_obj)   , INTENT(INOUT) :: texp,timp
      TYPE(bief_obj)   , INTENT(INOUT) :: t2_1,t2_2,t2_3,t3_1,t3_2
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     LOCAL VARIABLES
      INTEGER                     :: i
      DOUBLE PRECISION, PARAMETER :: eps=1.d-3
      DOUBLE PRECISION            :: cc
!
!-----------------------------------------------------------------------
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 0'
!
!     =======================================
!     PRELIMINARY COMPUTATIONS
!     =======================================
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 1'
!
!     BED SHEAR STRESS (TAUB-STOCKED IN T2_1==>2D TABLE)
!
      CALL taub_waqtel(cf,ro0,t2_1,npoin2,un,vn)
!
!     DEPOTION PROBABILITY (SED): STOCKED IN T2_2==>2D TABLE
!
      CALL depos_fx(t2_2,t2_1,tn%ADR(ind_ss)%P,taus,vitchu,npoin2)
!
!     EROSION FLUX (RS): STOCKED IN T2_3 ==> 2D TABLE
!
      CALL erosion_fx(t2_3,t2_1,tn%ADR(ind_sf)%P,taur,ero,1.d-10,
     &                npoin2)
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 4'
!
!
!     =======================================
!     LET'S NOW COMPUTE SOURCE TERMS
!     =======================================
!
!     FIRST TRACER: SUSPENDED LOAD [SS] (IND_SS)
!
!     BED SOURCES
!      DO I=1,NPOIN2
!        TEXP%ADR(IND_SS)%P%R(I)=T2_3%R(I)-T2_2%R(I)
!      ENDDO
!      CALL OVD('X=Y/Z   ',TEXP%ADR(IND_SS)%P%R,TEXP%ADR(IND_SS)%P%R,
!     &         ZPROP%R,0.D0,NPOIN2,2,0.D0,EPS                      )
      CALL os ('X=Y-Z   ', x=t2_1,y=t2_3,z=t2_2                     )
      CALL ovd('X=X+CY/Z',texp%ADR(ind_ss)%P%R,t2_1%R,
     &         zprop%R,1.d0,npoin2,2,0.d0,eps                       )
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 5'
!
!     SECOND TRACER: BED SEDIMENT [SF] (IND_SF)
!      warning: no advection neither diffusion for this tracer
!
      DO i=1,npoin2
!        TEXP%ADR(IND_SF)%P%R(I)=T2_2%R(I)-T2_3%R(I)
        texp%ADR(ind_sf)%P%R(i) = texp%ADR(ind_sf)%P%R(i) +
     &        t2_2%R(i) - t2_3%R(i)
      ENDDO
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 6'
!
!     THIRD TRACER: POLLUTANT DENSITY [C] (IND_C)
!
!     implicit part
      CALL os( 'X=X+C   ' ,x=timp%ADR(ind_c)%P,c=ccsedim            )
!     explicit part
!      CALL OS( 'X=CY    ' ,X=TEXP%ADR(IND_C)%P,Y=TN%ADR(IND_CSS)%P,
!     &                     C=KDESORP                                   )
      CALL os( 'X=CY    ' ,x=t3_2,y=tn%ADR(ind_css)%P     ,c=kdesorp )
      CALL os( 'X=X+Y   ' ,x=texp%ADR(ind_c)%P,y=t3_2                )
      cc =-kdesorp*cdistrib
!  warning: the following term causes divergence of the code, it should
!           be traited implicitly- it is commented: to be investigated
!           more in depth
!      CALL OS( 'X=X+CYZ ' ,X=TEXP%ADR(IND_C)%P,Y=TN%ADR(IND_C)%P,
!     &                     Z=TN%ADR(IND_SS)%P             ,C=CC       )
!
!
!     FORTH TRACER: ABSORBED POLLUTANT BY SUSPENDED LOAD [CSS] (IND_CSS)
!
!     implicit part
      CALL os( 'X=X+C   ' ,x=timp%ADR(ind_css)%P,c=ccsedim          )
!     explicit part
!      CALL OS( 'X=-Y    ' ,X=TEXP%ADR(IND_CSS)%P,Y=TEXP%ADR(IND_C)%P )
      CALL os( 'X=X-Y   ' ,x=texp%ADR(ind_css)%P,y=t3_2 )
      DO i=1,npoin2
        t3_1%R(i) = t2_3%R(i)*tn%ADR(ind_csf)%P%R(i) -
     &              t2_2%R(i)*tn%ADR(ind_css)%P%R(i)
      ENDDO
      CALL ovd('X=Y/Z   ' ,t3_1%R,t3_1%R,
     &                     zprop%R,0.d0,npoin2,2,0.d0,eps           )
      CALL os( 'X=X+Y   ' ,x=texp%ADR(ind_css)%P,y=t3_1             )
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 8'
!
!     FIFTH TRACER: ABSORBED POLLUTANT BY BED SEDIMENT [CFF] (IND_CSF)
!
      DO i=1,npoin2
        texp%ADR(ind_csf)%P%R(i) = texp%ADR(ind_csf)%P%R(i) +
     &                           t2_2%R(i)*tn%ADR(ind_css)%P%R(i) -
     &                           t2_3%R(i)*tn%ADR(ind_csf)%P%R(i)
      ENDDO
      CALL os( 'X=X+CY  ' ,x=texp%ADR(ind_csf)%P,y=tn%ADR(ind_csf)%P,
     &                     c=-ccsedim                               )
!
      IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 9'
!
!     MASS BALANCE: MASS ADDED BY EXPLICIT TERMS
!                   (IMPLICIT PART IS ADDED IN CVDFTR)
!
!-----------------------------------------------------------------------
!
      RETURN
      END