calcs2d_micropol.f

Go to the documentation of this file.

!                     ***************************
                      SUBROUTINE calcs2d_micropol
!                     ***************************
!
     & (npoin,tn,texp,timp,hprop,cf,un,vn,
     &  t1,t2,t3,t4)
!
!***********************************************************************
! WAQTEL   V8P1
!***********************************************************************
!
!brief    COMPUTES SOURCE TERMS FOR MICROPOL WAQ PROCESS
!          WAQ PROCESS OF CODE_TRACER (MASCARET SYSTEM)
!
!history  R. ATA
!+        21/09/2014
!+        V7P0
!+       CREATION (VOID)
!history  R. ATA
!+        28/09/2015
!+        V7P1
!+       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
!
!history  C.-T. PHAM
!+        20/11/2019
!+        V8P1
!+        The calculation of the 4th ad 5th equations is with RS/SF
!+        and SEDP not RS and SED=SEDP*SS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CCSEDIM        |-->| CONSTANT OF EXPONENTIAL DESINTEGRATION
!| CDISTRIB       |-->| COEFFICIENT OF DISTRIBUTION (KD)
!| ERO            |-->| EROSION RATE
!| KDESORP        |-->| KINETIC CONSTANT OF  DESORPTION
!| NPOIN          |-->| TOTAL NUMBER OF MESH NODES
!| NTRAC          |-->| NUMBER OF TRACERS
!| 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_calcs2d_micropol => calcs2d_micropol
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      INTEGER          , INTENT(IN   ) :: npoin
      TYPE(bief_obj)   , INTENT(IN   ) :: tn,hprop,cf,un,vn
      TYPE(bief_obj)   , INTENT(INOUT) :: texp,t1,t2,t3,t4,timp
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     LOCAL VARIABLES
      DOUBLE PRECISION, PARAMETER :: eps=1.d-3
      DOUBLE PRECISION            :: cc
!
!-----------------------------------------------------------------------
!
!
!     =======================================
!     PRELIMINARY COMPUTATIONS
!     =======================================
!
!     BED SHEAR STRESS (TAUB-STORED IN T1)
!
      CALL taub_waqtel(cf,ro0,t1,npoin,un,vn)
!
!     DEPOSITION PROBABILITY (SED): STORED IN T2, NOT SEDP
!
      CALL depos_fx(t2,t1,tn%ADR(ind_ss)%P,taus,vitchu,npoin)
!
!     SEDP = SED/SS: STORED IN T4, NOT SED
!
      CALL os ('X=C     ', x=t3,c=1.d0)
      CALL depos_fx(t4,t1,t3,taus,vitchu,npoin)
!
!     EROSION FLUX (RS): STORED IN T3
!
      CALL erosion_fx(t3,t1,tn%ADR(ind_sf)%P,taur,ero,1.d-10,npoin)
!
!
!     =======================================
!     LET'S NOW COMPUTE SOURCE TERMS
!     =======================================
!
!     FIRST TRACER: SUSPENDED LOAD [SS] (IND_SS)
!                                               RS   SED
      CALL os ('X=Y-Z   ', x=t1,              y=t3,z=t2              )
      CALL ovd('X=X+CY/Z', texp%ADR(ind_ss)%P%R,t1%R,hprop%R,
     &                     1.d0, npoin,2,0.d0,eps  )
!
!     SECOND TRACER: BED SEDIMENT [SF] (IND_SF)
!      warning: no advection neither diffusion for this tracer
!                                SED   RS
      CALL os ('X=Y-Z   ', x=t1,y=t2,z=t3              )
      CALL os ('X=X+Y   ', x=texp%ADR(ind_sf)%P, y=t1  )
!
!     THIRD TRACER: POLLUTANT DENSITY [C] (IND_C)
!
!     implicit part
!     CCSEDIM = CONSTANT L IN THE DOC      
!      CALL OS( 'X=X+C   ' ,X=TIMP%ADR(IND_C)%P,C=CCSEDIM             )
      CALL os( 'X=X+CY  ' ,x=timp%ADR(ind_c)%P,y=hprop,c=-ccsedim)
!     explicit part
      CALL os( 'X=CY    ' ,x=t1,y=tn%ADR(ind_css)%P,c=kdesorp        )
      cc =-kdesorp*cdistrib
      CALL os( 'X=X+CYZ ' ,x=t1,y=tn%ADR(ind_c)%P,
     &                     z=tn%ADR(ind_ss)%P             ,c=cc       )
      CALL os( 'X=X+Y   ' ,x=texp%ADR(ind_c)%P,y=t1                  )
!
!     FOURTH TRACER: ADSORBED POLLUTANT BY SUSPENDED LOAD [CSS] (IND_CSS)
!
!     implicit part
!      CALL OS( 'X=X+C   ' ,X=TIMP%ADR(IND_CSS)%P,C=CCSEDIM             )
      CALL os( 'X=X+CY  ' ,x=timp%ADR(ind_css)%P,y=hprop,c=-ccsedim)
!     explicit part
      CALL os( 'X=X-Y   ' ,x=texp%ADR(ind_css)%P,y=t1                  )
!      CALL OS( 'X=YZ    ' ,X=T4,Y=T3,Z=TN%ADR(IND_CSF)%P               )
!      CALL OS( 'X=X+CYZ ' ,X=T4,Y=TN%ADR(IND_CSS)%P,Z=T2,C=-1.D0       )
!      CALL OVD('X=Y/Z   ' ,T4%R,T4%R,
!     &                     HPROP%R,0.D0,NPOIN,2,0.D0,EPS)
!      CALL OS( 'X=X+Y   ' ,X=TEXP%ADR(IND_CSS)%P,Y=T4                  )
!                                 RS
      CALL os( 'X=Y     ' ,x=t1,y=t3)
      CALL ovd('X=CXY/Z ' ,t1%R,tn%ADR(ind_csf)%P%R,tn%ADR(ind_sf)%P%R,
     &                     1.d0, npoin,2,0.d0,eps)
!                                 SEDP
      CALL os( 'X=X-YZ  ' ,x=t1,y=t4,z=tn%ADR(ind_css)%P)
      CALL ovd('X=X+CY/Z', texp%ADR(ind_css)%P%R,t1%R,hprop%R,
     &                     1.d0, npoin,2,0.d0,eps  )

!
!     FIFTH TRACER: ADSORBED POLLUTANT BY BED SEDIMENT [CFF] (IND_CSF)
!
!      CALL OS( 'X=X+YZ  ' ,X=TEXP%ADR(IND_CSF)%P,Y=TN%ADR(IND_CSS)%P,
!     &                     Z=T2                                        )
!      CALL OS( 'X=X+CYZ ' ,X=TEXP%ADR(IND_CSF)%P,Y=TN%ADR(IND_CSF)%P,
!     &                     Z=T3,C=-1.D0                                )
!      CALL OS( 'X=X+CY  ' ,X=TEXP%ADR(IND_CSF)%P,Y=TN%ADR(IND_CSF)%P,
!     &                     C=-CCSEDIM                                  )
!     implicit part
      CALL os( 'X=X+CY  ' ,x=timp%ADR(ind_csf)%P,y=hprop,c=-ccsedim)
!     explicit part
      CALL os( 'X=X-Y   ' ,x=texp%ADR(ind_csf)%P,y=t1)
!
!-----------------------------------------------------------------------
!
      RETURN
      END