reaer.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE reaer
!                   ****************
!
     &(formk2,k2,k22,npoin2,nplan,un,vn,h,eps)
!
!***********************************************************************
! WAQTEL   V8P1
!***********************************************************************
!
!brief    COMPUTES THE COEFFICIENT OF REAERATION K2
!
!
!
!history  R. ATA (LNHE)
!+        02/09/2015
!+        V7P1
!+
!history  R. ATA (LNHE)
!+        22/02/2016
!+        V7P2
!+  ADD CONVERSION SEC TO DAY
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| EPS            |-->| TO AVOID DIVISION BY 0
!| FORMK2         |-->| WHICH FROMULA TO COMPUTE K2
!| H              |-->| WATER DEPTH
!| K2             |<--| COEFFICIENT OF REAERATION
!| K22            |-->| COEFFICIENT OF REAERATION IF IT IS CONSTANT
!| NPLAN          |-->| NUMBER OF VERTICAL LAYERS
!| NPOIN2         |-->| NUMBER OF MESH NODES OF THE SURFACE LAYER
!| UN,VN          |-->| VELOCITY COMPONENTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_special
      USE declarations_waqtel,ONLY: sectoday
      USE interface_waqtel, ex_reaer => reaer
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          ,INTENT(IN   ) :: FORMK2,NPOIN2,NPLAN
      DOUBLE PRECISION ,INTENT(IN   ) :: EPS,K22
      TYPE(bief_obj)   ,INTENT(IN   ) :: H,UN,VN
      TYPE(bief_obj)   ,INTENT(INOUT) :: K2
      INTRINSIC abs,sqrt,max
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
!     LOCAL VARIABLES
      INTEGER           I,J
      DOUBLE PRECISION  UNORM,PJ
!
!
      IF(formk2.EQ.0)THEN ! ==> CONSTANT K2
        CALL os('X=C     ',x=k2,c=k22)
      ELSE  ! ==> VARIABLE K2
        DO j =1,npoin2
!
          i=(nplan-1)*npoin2+j
          unorm = sqrt(un%R(i)**2+vn%R(i)**2) !  GENERALIZATION OF U IN 1D
!         HEAD OR WATER DEPTH TO REPLACE RH ?
          pj = h%R(i)  !ZF%R(I)+H%R(I)+UNORM**2/DEUXG
!         FORMULA OF THE TENESSEE VALLEY AUTHORITY
          IF( formk2.EQ.1 ) THEN
            k2%R(j) = 5.23d0*unorm* max(h%R(i),eps)**(-1.67d0)
!         FORMULA OF OWENS ET AL.
          ELSEIF(formk2.EQ.2)THEN
            k2%R(j) = 5.33d0*(unorm**0.67d0)*
     &              max(h%R(i),eps)**(-1.85d0)
!         FORMULA OF CHURCHILL ET AL.
          ELSEIF(formk2.EQ.3)THEN
            k2%R(j) = 0.746d0 * (unorm**2.695d0) /
     &      (max(h%R(i),eps)**3.085d0 *
     &       max(abs(pj),eps)**0.823d0) ! VERIFY THE FORMULA, SOME DOUBT ?
!         FORMULA OF O CONNOR & DOBBINS'
          ELSEIF(formk2.EQ.4)THEN
            k2%R(j) = (3.90d0 * unorm**0.5d0 ) /
     &               max(h%R(i),eps)**1.5d0
!         FORMULA OF ?? INVISIBLE MAN :) : IT SEEMS TO BE A COMBINATION OF THE 3 LAST FORMULA ?!
          ELSEIF( formk2.EQ.5 ) THEN
            IF( h%R(i).LE.0.6d0 ) THEN
              k2%R(j) = 5.33d0 * (unorm**0.67d0) *
     &                max(h%R(i),eps)**(-1.85d0)
            ELSEIF (h%R(i).LT.(12.d0*unorm-6.6d0)) THEN
              k2%R(j) =  0.746d0*(unorm**2.695d0)/
     &                (max(h%R(i), eps)**3.085d0 *
     &                 max(abs(pj),eps)**0.823d0)
            ELSE
              k2%R(j) = 3.90d0 * (unorm**0.5d0)/
     &                max(h%R(i),eps)**1.5d0
            ENDIF
          ELSE
            WRITE(lu,111)formk2
            CALL plante(1)
            stop
          ENDIF
        ENDDO
      ENDIF
!
!     CONVERT DAY TO SEC
!
!      CALL OS('X=CX    ',X=K2,C=SECTODAY)
!
!     ERROR MESSAGES
!
!
111   FORMAT(1x,'REAR.F: K2 FORMULA :',i3,/,1x, 'NOT AVAILABLE')
!
      RETURN
      END