positive_depths.f

Go to the documentation of this file.

!                   **************************
                    SUBROUTINE positive_depths
!                   **************************
!
     &(t1,t2,t3,t4,h,hn,mesh,flodel,compute_flodel,flbor,dt,
     & unsv2d,npoin,gloseg1,gloseg2,nbor,nptfr,
     & smh,yasmh,pluie,rain,optsou,flulim,limpro,hbor,kdir,info,
     & flopoint,namecode,option,nitmax,doflulim,flulimebe,doflulimebe)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    SUPPRESSES NEGATIVE DEPTHS BY A LIMITATION OF FLUXES.
!
!history  J-M HERVOUET (LNHE)
!+        12/03/2010
!+        V6P0
!+   ADDED COMPUTE_FLODEL
!
!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 (LNHE)
!+        16/07/2012
!+        V6P2
!+   The value of transmitted fluxes is given back in FLODEL in all
!+   cases, YAFLODEL no longer used (but left for compatibility)
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        30/05/2013
!+        V6P2
!+   Argument NITMAX added.
!
!history  J-M HERVOUET (LNHE)
!+        30/12/2013
!+        V7P0
!+   Argument YAFLODEL removed.
!
!history  J-M HERVOUET (LNHE)
!+        18/09/2014
!+        V7P0
!+   Now positive sources are first added, then the transfers between
!+   points are done, then the negative sources are treated. In this way
!+   the transfers cannot be corrupted by an initial negative depth.
!+   + one NCSIZE.GT.0 changed into NCSISE.GT.1.
!
!history  J-M HERVOUET (LNHE)
!+        19/06/2015
!+        V7P1
!+   Adaptation to new coefficient IFAC instead of FAC.
!+   Option 1 removed. FLULIM now equal on either side of a segment
!+   in parallel. FLODEL differently shared at the exit.
!
!history  J-M HERVOUET (LNHE)
!+        25/03/2016
!+        V7P2
!+   OPTION has a new meaning (2 was the only previous possibility)
!+   Now : 1 : positive depths ensured with an EBE approach (new !)
!+             that is necessary for the ERIA advection scheme.
!+         2 : positive depths ensured with an edge-based approach (old)
!+   Option 1 will be mandatory with the new advection scheme 15.
!
!history  J-M HERVOUET (LNHE)
!+        16/07/2016
!+        V7P2
!+   In ERIA scheme, now forbidding that two fluxes on either side of
!+   a segment be of opposite sign (one is cancelled and the other is
!+   reduced). The final flux limitation is thus in the range [0,1] and
!+   OPTION=1 (i.e. TREATMENT OF NEGATIVE DEPTHS=3) is now also possible
!+   for the LIPS scheme.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        09/09/2016
!+        V7P2
!+   Swapping rain and boundary fluxes at the end of the algorithm. It
!+   caused a wrong correction of FLBOR in case of evaporation.
!+   Adding rain (it was mixed before with sources).
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        15/09/2016
!+        V7P2
!+   Limiting evaporation to avoid negative depths.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        06/11/2016
!+        V7P3
!+   This subroutine can now provide a flux limitation in the form of
!+   FLULIM(NSEG) or FLULIMEBE(NELMAX,3). This depends on OPTION
!+   However, FLULIMEBE can be produced only with OPTION 1.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        19/11/2016
!+        V7P3
!+   Subroutine split into POSITIVE_DEPTHS_ERIA and POSITIVE_DEPTHS_NERD
!+   The remaining part is just the bifurcation.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| COMPUTE_FLODEL |-->| IF YES, COMPUTE FLODEL HERE
!| DOFLULIM       |-->| OPTIONAL, IF YES DOES ARRAY FLULIM
!| DOFLULIMEBE    |-->| OPTIONAL, IF YES DOES ARRAY FLULIMEBE
!| DT             |-->| TIME STEP
!| FLBOR          |<->| BOUNDARY FLUXES
!| FLODEL         |<->| FLUXES GIVEN BY SEGMENT
!|                |   | MAY BE COMPUTED HERE (SEE COMPUTE-FLODEL)
!|                |   | OR SIMPLY GIVEN. AT THE EXIT, THE REAL FLUX
!|                |   | TRANSMITTED IS GIVEN BACK.
!| FLOPOINT       |-->| FLUXES GIVEN BY POINTS (ELEMENT BY ELEMENT)
!| FLULIM         |<--| PER SEGMENT: PERCENTAGE OF FLUX THAT HAS NOT
!|                |   | BEEN TRANSMITTED AT THE END OF THE ALGORITHM
!| FLULIMEBE      |<--| AS FLULIM BUT GIVEN PER ELEMENT
!|                |   | HENCE FLULIMEBE(NELMAX,3)
!| GLOSEG1        |-->| FIRST POINT OF SEGMENTS
!| GLOSEG2        |-->| SECOND POINT OF SEGMENTS
!| H              |<->| NEW DEPTH
!| HBOR           |-->| PRESCRIBED DEPTHS AT BOUNDARIES
!| HN             |-->| OLD DEPTH
!| INFO           |-->| IF YES, PRINTING INFORMATION ON LISTING
!| KDIR           |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITION
!| LIMPRO         |-->| TYPE OF BOUNDARY CONDITIONS
!|                |   | IF EQUAL TO KDIR: PRESCRIBED DEPTH.
!| MESH           |<->| MESH STRUCTURE
!| NAMECODE       |-->| NAME OF CALLING CODE (SISYPHE, TELEMEC2D, ETC.)
!| NBOR           |-->| GLOBAL NUMBERS OF BOUNDARY POINTS
!| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
!| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| OPTION         |-->| OPTION OF ALGORITHM FOR EDGE-BASED ADVECTION
!|                |   | 1: TRIANGLES, COMPATIBLE WITH ERIA
!|                |   | 2: SEGMENTS, COMPATIBLE WITH NERD
!| OPTSOU         |-->| OPTION FOR SOURCES 1: NORMAL 2: DIRAC
!| PLUIE          |-->| RAIN IN A BIEF_OBJ, IN M/S.
!| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
!| SMH            |-->| SOURCE TERMS
!| T1             |-->| WORK ARRAY
!| T2             |-->| WORK ARRAY
!| T3             |-->| WORK ARRAY
!| T4             |-->| WORK ARRAY
!| UNSV2D         |-->| INVERSE OF INTEGRAL OF BASIS FUNCTIONS
!| YASMH          |-->| IF(YES) SMH MUST BE TAKEN INTO ACCOUNT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_positive_depths => positive_depths
      USE interface_parallel
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NPOIN,NPTFR,OPTSOU,KDIR,OPTION
      INTEGER, INTENT(IN)             :: NITMAX
      INTEGER, INTENT(IN)             :: GLOSEG1(*),GLOSEG2(*)
      INTEGER, INTENT(IN)             :: NBOR(nptfr)
      INTEGER, INTENT(IN)             :: LIMPRO(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,HBOR(nptfr)
      DOUBLE PRECISION, INTENT(INOUT) :: FLULIM(*)
      TYPE(bief_mesh),INTENT(INOUT)   :: MESH
      DOUBLE PRECISION, INTENT(INOUT) :: FLOPOINT(mesh%nelmax,3)
      TYPE(bief_obj), INTENT(INOUT)   :: T1,T2,T3,T4,FLODEL,H,FLBOR
      TYPE(bief_obj), INTENT(INOUT)   :: PLUIE
      TYPE(bief_obj), INTENT(IN)      :: UNSV2D,HN,SMH
      LOGICAL, INTENT(IN)             :: YASMH,INFO,RAIN,COMPUTE_FLODEL
      CHARACTER(LEN=24)               :: NAMECODE
      LOGICAL, INTENT(IN), OPTIONAL   :: DOFLULIM,DOFLULIMEBE
      DOUBLE PRECISION, INTENT(INOUT),
     &                      OPTIONAL, TARGET :: flulimebe(mesh%NELMAX,3)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL MAKEFLULIM,MAKEFLULIMEBE
      DOUBLE PRECISION, POINTER :: PT_FLULIMEBE(:,:)
!
!-----------------------------------------------------------------------
!
      IF(PRESENT(doflulim)) THEN
        makeflulim=doflulim
      ELSE
!       DEFAULT VALUE (FLULIM WAS ALWAYS DONE BEFORE INTRODUCING DOFLULIM)
        makeflulim=.true.
        ENDIF
      IF(PRESENT(doflulimebe)) THEN
        makeflulimebe=doflulimebe
      ELSE
!       DEFAULT VALUE (FLULIMEBE WAS NOT DONE BEFORE INTRODUCING DOFLULIMEBE)
        makeflulimebe=.false.
      ENDIF
!
      IF(makeflulimebe.AND..NOT.PRESENT(flulimebe)) THEN
        WRITE(lu,*) 'POSITIVE_DEPTHS: '
        WRITE(lu,*) 'MISSING ARGUMENT FLULIMEBE'
        CALL plante(1)
        stop
      ELSEIF(PRESENT(flulimebe)) THEN
        pt_flulimebe=>flulimebe
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(option.EQ.1) THEN
        CALL positive_depths_eria(t1,t2,t3,t4,h,hn,mesh,flodel,
     &                            compute_flodel,flbor,dt,
     &                            unsv2d,npoin,gloseg1,gloseg2,
     &                            nbor,nptfr,smh,yasmh,pluie,rain,
     &                            optsou,limpro,hbor,kdir,info,
     &                            flopoint,namecode,nitmax,
     &                            makeflulimebe,pt_flulimebe)
      ELSEIF(option.EQ.2) THEN
        CALL positive_depths_nerd(t1,t2,t4,h,hn,mesh,flodel,
     &                            compute_flodel,flbor,dt,
     &                            unsv2d,npoin,gloseg1,gloseg2,
     &                            nbor,nptfr,smh,yasmh,pluie,rain,
     &                            optsou,flulim,limpro,hbor,kdir,info,
     &                            flopoint,namecode,nitmax,makeflulim)
      ELSE
        WRITE(lu,*) 'UNKNOWN OPTION IN POSITIVE_DEPTHS: ',option
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END