bilant.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE bilant
!                   *****************
!
     &(h,work2,work3,dt,lt,nit,info,
     & t,agglot,massou,mastr0,mastr2,masten,
     & mastou,msk,maskel,mesh,
     & numliq,nfrliq,nptfr,nametrac,flbortra,mass_rain,train,
     & mastrain)
!
!***********************************************************************
! BIEF   V7P0
!***********************************************************************
!
!brief    COMPUTES THE MASS BALANCE FOR THE TRACER.
!
!history  J-M HERVOUET (LNHE)     ; C MOULIN (LNH)
!+        10/06/08
!+        V5P9
!+
!
!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  J-M HERVOUET (EDF LAB, LNHE)
!+        30/09/20148
!+        V7P0
!+   Some ABS put in the relative accuracy formulas because the mass of
!+   a tracer can be negative (case of vorticity).
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AGGLOT         |-->| MASS-LUMPING ON TRACER
!| DT             |-->| TIME-STEP
!| FLBORTRA       |-->| TRACER FLUXES AT BOUNDARIES
!| H              |-->| DEPTH AT TIME N+1.
!| INFO           |-->| LOGICAL, IF YES, PRINTING INFORMATION ON LISTING
!| LT,NIT         |-->| TIME STEP NUMBER, TOTAL NUMBER OF STEPS.
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MASS_RAIN      |<--| MASS OF WATER ADDED BY RAIN OR EVAPORATION
!| MASSOU         |-->| MASS OF TRACER BROUGTH BY SOURCE TERM
!| MASTEN         |<--| WATER MASS ENTERED THROUGH BOUNDARIES
!| MASTRAIN       |<->| TOTAL WATER MASS ENTERED THROUGH BOUNDARIES
!| MASTOU         |<--| WATER MASS CREATED BY SOURCE TERMS
!| MASTR0         |<--| INITIAL TRACER MASS
!| MASTR2         |<--| CURRENT TRACER MASS
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| NAMETRAC       |-->| NAMES OF TRACERS
!| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NUMLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
!| T              |-->| TRACER AT TIME T(N+1)
!| TRAIN          |-->| VALUE OF TRACER IN THE RAIN
!| WORK2          |<->| WORK ARRAY
!| WORK3          |<->| WORK ARRAY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_bilant => bilant
!
      USE declarations_special
      USE interface_parallel, ONLY : p_sum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)            :: LT,NIT,NFRLIQ,NPTFR
      INTEGER, INTENT(IN)            :: NUMLIQ(nptfr)
      DOUBLE PRECISION, INTENT(IN)   :: DT,MASSOU,AGGLOT,MASS_RAIN,TRAIN
      DOUBLE PRECISION, INTENT(INOUT):: MASTRAIN
      LOGICAL, INTENT(IN)            :: INFO,MSK
      TYPE(bief_obj), INTENT(INOUT)  :: WORK2,WORK3
      TYPE(bief_obj), INTENT(IN)     :: H,T,MASKEL
      TYPE(bief_obj), INTENT(IN)     :: FLBORTRA
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
      DOUBLE PRECISION, INTENT(INOUT):: MASTR0,MASTR2,MASTEN,MASTOU
      CHARACTER(LEN=32), INTENT(IN)  :: NAMETRAC
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IFRLIQ,IELMT,IELMH
!
      DOUBLE PRECISION ERREUT,PERDUE,FLUXT,MASBOR,RELATI,DENOM,MASTR1
      DOUBLE PRECISION MASRAI
!     300 IS HERE MAXFRO, THE MAXIMUM NUMBER OF LIQUID BOUNDARIES
      DOUBLE PRECISION FLT_BOUND(300)
!
      INTRINSIC abs,max
!
!-----------------------------------------------------------------------
!
      ielmt = t%ELM
      ielmh = h%ELM
!
!-----------------------------------------------------------------------
!
!  COMPATIBLE COMPUTATION OF THE TRACER QUANTITY AT TIME N+1:
!  TAKES MASS-LUMPING INTO ACCOUNT BUT REQUIRES AGGLOC=AGGLOT
!
      IF(lt.NE.0) mastr1 = mastr2
!
      CALL vector(work2,'=','MASVEC          ',ielmt,
     &            1.d0-agglot,t,t,t,t,t,t,mesh,msk,maskel)
!     H IS GIVEN HERE FOR A DUMMY STRUCTURE
      CALL vector(work3,'=','MASBAS          ',ielmt,
     &                 agglot,h,h,h,h,h,h,mesh,msk,maskel)
!
      CALL os('X=X+YZ  ',x=work2,y=work3,z=t)
!
      mastr2 = dots(work2,h)
      IF(ncsize.GT.1) mastr2=p_sum(mastr2)
!
      IF(lt.EQ.0) THEN
        mastr0   = mastr2
        mastr1   = mastr2
        masten   = 0.d0
        mastou   = 0.d0
        mastrain = 0.d0
      ENDIF
!
!=======================================================================
!
!   COMPUTES THE FLUXES (MISSES THE DIFFUSION FLUX,... INVESTIGATE)
!
!=======================================================================
!
      fluxt=0.d0
!
      IF(lt.GT.0.AND.nfrliq.GT.0) THEN
        DO ifrliq=1,nfrliq
          flt_bound(ifrliq)=0.d0
        ENDDO
        IF(nptfr.GT.0) THEN
          DO i=1,nptfr
!           NOTE: COULD DEFINE FLUX_BOUNDARIES BETWEEN 0 AND NFRLIQ
            ifrliq=numliq(i)
            IF(ifrliq.GT.0) THEN
!             FLBORTRA MUST NOT BE ASSEMBLED IN PARALLEL MODE
              flt_bound(ifrliq)=flt_bound(ifrliq)+flbortra%R(i)
            ENDIF
          ENDDO
        ENDIF
        IF(ncsize.GT.1) THEN
          DO ifrliq=1,nfrliq
            flt_bound(ifrliq)=p_sum(flt_bound(ifrliq))
          ENDDO
        ENDIF
        DO ifrliq=1,nfrliq
          fluxt=fluxt+flt_bound(ifrliq)
        ENDDO
      ENDIF
!
!=======================================================================
!
!     COMPUTES THE FLUXES AT THE LIQUID BOUNDARIES
!
      masten = masten - fluxt * dt
      mastou = mastou + massou
!
!=======================================================================
!
!     COMPUTES THE TRACER FLUXES THRU THE WALLS (FLUX LAW)
!
!     TEMPORARY, TO BE CODED UP
      masbor = 0.d0
!
!=======================================================================
!
!     COMPUTES THE TRACER MASS BROUGHT BY THE RAIN
!     MAX(...,0.D0) BECAUSE EVAPORATION IS NOT TAKEN INTO ACCOUNT
!
!     WILL WORK ONLY IF TRAIN AND RAIN CONSTANT ON ALL THE DOMAIN...
      masrai = max(mass_rain,0.d0) * train
      IF(ncsize.GT.1) masrai=p_sum(masrai)
      mastrain = mastrain + masrai
!
!=======================================================================
!
!     COMPUTES THE ERROR ON THE MASS FOR THIS TIMESTEP
!
      erreut = mastr1 + massou + masrai - mastr2 - dt*fluxt
!
!=======================================================================
!
!     PRINTOUTS :
!
      IF(info) THEN
!
!-----------------------------------------------------------------------
!
!     PRINTOUTS FOR THE TRACER
!
        WRITE(lu,*)
        WRITE(lu,*) '                           BALANCE OF ',
     &    trim(nametrac(1:16)),' (UNIT: ',trim(nametrac(17:32)),' * M3)'
!
        IF(lt.EQ.0) THEN
          WRITE(lu,2090) mastr0
        ELSE
          WRITE(lu,2160) mastr1,mastr2
          IF(nfrliq.GT.0) THEN
            DO ifrliq=1,nfrliq
              WRITE(lu,2110) ifrliq,-flt_bound(ifrliq)
            ENDDO
          ENDIF
          IF(abs(massou).GT.1.d-8) THEN
            WRITE(lu,2113) massou
          ENDIF
          IF(abs(masrai).GT.1.d-8) THEN
            WRITE(lu,2166) masrai
          ENDIF
          WRITE(lu,2165) erreut
!         ABS BECAUSE THE MASS OF A TRACER CAN BE NEGATIVE
!         EXAMPLE : VORTICITY
          denom = max(abs(mastr1),abs(mastr2),abs(fluxt*dt),
     &                abs(masrai),abs(massou))
          IF(denom.GT.1.d-8) THEN
            erreut = erreut / denom
            WRITE(lu,2120) erreut
          ENDIF
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  FINAL MASS BALANCE
!
      IF(lt.EQ.nit) THEN
!
        WRITE(lu,*)
        WRITE(lu,*) '                     FINAL BALANCE OF ',
     &    trim(nametrac(1:16)),' (UNIT: ',trim(nametrac(17:32)),' * M3)'
!
        perdue = mastr0+masten+masbor+mastou+mastrain-mastr2
        denom = max(abs(mastr0),abs(mastr2),abs(masten),
     &              abs(mastou),abs(mastrain))
        WRITE(lu,2160) mastr0,mastr2
        IF(abs(masten).GT.1.d-8) WRITE(lu,2161) masten
        IF(abs(mastou).GT.1.d-8) WRITE(lu,2164) mastou
        IF(abs(mastrain).GT.1.d-8) WRITE(lu,2166) mastrain
        WRITE(lu,2165) perdue
        IF(denom.GT.1.d-8) THEN
          relati = perdue / denom
          WRITE(lu,2120) relati
        ENDIF
        WRITE(lu,*)
!
      ENDIF
!
!  END OF THE PRINTOUTS :
!
!=======================================================================
!
!  FORMATS :
!
2090  FORMAT(  5x,'INITIAL QUANTITY OF TRACER   :',g16.7)
2110  FORMAT(  5x,'BOUNDARY ',1i3,' FLUX:         ',g16.7,
     &          ' ( >0 : ENTERING  <0 : EXITING )')
2113  FORMAT(  5x,'QUANTITY CREATED BY SOURCE TERM   :   ' , g16.7 )
2120  FORMAT(  5x,'RELATIVE ERROR                    : ',g16.7)
2160  FORMAT(/,5x,'INITIAL QUANTITY                  : ',g16.7,
     &       /,5x,'FINAL QUANTITY                    : ',g16.7)
2161  FORMAT(  5x,'QUANTITY ENTERED THROUGH LIQ. BND.: ',g16.7,
     &            '  ( IF <0 EXIT )')
2164  FORMAT(  5x,'QUANTITY CREATED BY SOURCE TERM   : ',g16.7)
2166  FORMAT(  5x,'MASS BROUGHT BY THE RAIN          : ',g16.7)
2165  FORMAT(  5x,'TOTAL QUANTITY LOST               : ',g16.7)
!
!=======================================================================
!
      RETURN
      END