calcs3d_o2.f

Go to the documentation of this file.

!                     *********************
                      SUBROUTINE calcs3d_o2
!                     *********************
!
     & (npoin3,npoin2,nplan,wattemp,o2satu,demben,formk2,k1,k44,k22,
     &  photo,resp,tn,texp,timp,t31,t32,t21,hprop,zprop,un,vn)
!
!***********************************************************************
! WAQTEL   V8P1
!***********************************************************************
!
!brief    COMPUTES SOURCE TERMS FOR O2 WAQ PROCESS FOR 3D CASE
!
!history  R. ATA
!+        21/02/2016
!+        V7P2
!+       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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DEMBEN        |-->| BENTIC DEMAND
!| DT            |-->| TIME STEP
!| FORMK2        |-->| FORMULA FOR COMPUTINK K2
!| HPROP         |-->| PROPAGATION DEPTH
!| K1            |-->| CONST. OF DEGRADATION KINETIC OF ORGANIC LOAD
!| K22           |<--| CONST. OF REAERATION
!| K44           |-->| CONST. OF NITRIFICATION KINETIC
!| MASSOU        |<--| MASS OF TRACER ADDED BY SOURCE TERM
!| NPOIN         |-->| NUMBER OF NODES IN THE MESH
!| O2SATU        |<--| O2 SATURATION DENSITY OF WATER (CS)
!| PHOTO         |-->| PHOTOSYNTHESIS
!| WATTEMP       |-->| TEMPERATURE OF WATER (DEG CELSIUS)
!| T1,T2         |-->| WORKING ARRAYS
!| TEXP          |<--| EXPLICIT SOURCE TERM.
!| TIMP          |<--| IMPLICIT SOURCE TERM.
!| TN            |-->| TRACERS AT TIME N
!| UN,VN         |-->| VELOCITY COMPONENTS AT TIME N
!| ZPROP         |-->| PROPAGATION ELEVATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

      USE bief
      USE declarations_special
      USE declarations_waqtel,ONLY: formcs,k2,sectoday,
     &  ind_t,ind_o2,ind_ol,ind_nh4
      USE interface_parallel
      USE interface_waqtel, ex_calcs3d_o2 => calcs3d_o2
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      INTEGER          , INTENT(IN   ) :: npoin2,npoin3,nplan
      INTEGER          , INTENT(IN   ) :: formk2
      DOUBLE PRECISION , INTENT(IN   ) :: demben,wattemp
      DOUBLE PRECISION , INTENT(IN   ) :: photo,resp,k1,k44
      DOUBLE PRECISION , INTENT(INOUT) :: o2satu,k22
      TYPE(bief_obj)   , INTENT(IN   ) :: tn,hprop,un,vn,zprop
      TYPE(bief_obj)   , INTENT(INOUT) :: texp,timp,t31,t32,t21
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     LOCAL VARIABLES
      INTEGER                     :: i,j
      DOUBLE PRECISION, PARAMETER :: eps=1.d-3
      DOUBLE PRECISION, PARAMETER :: corr1=1.065d0
      DOUBLE PRECISION, PARAMETER :: corr2=1.0241d0
      DOUBLE PRECISION            :: bencorr,pmoinr,conv_k1,conv_k44
      DOUBLE PRECISION            :: power
      INTRINSIC max
!
      pmoinr  = (photo-resp)*sectoday
!     CONVERT DAYS TO SECONDS
      conv_k44= k44   *sectoday
      conv_k1 = k1    *sectoday
      bencorr = demben*sectoday
      power   = wattemp-20.d0
      DO i=1,npoin3
        IF( ind_t.GT.0 ) power=tn%ADR(ind_t)%P%R(i)-20.d0
!       CORR2T AND BENCOR STOCKED HERE IN T31,T32
        t31%R(i)=corr2**power
        t32%R(i)=bencorr*(corr1**power)
      ENDDO
!
!     COMPUTE CS (O2SATU, STOCKED IN T21)
!
      IF( ind_t.EQ.0 )THEN
        CALL satur_o2(o2satu,formcs,wattemp,eps)
        CALL os('X=C     ',x=t21,c=o2satu       )
      ELSE
        DO i=1,npoin2
          j=(nplan-1)*npoin2+i
          CALL satur_o2(t21%R(i),formcs,tn%ADR(ind_t)%P%R(j),eps)
        ENDDO
      ENDIF
!
!     COMPUTE K2
!
      CALL reaer(formk2,k2,k22,npoin2,nplan,un,vn,hprop,eps)
!     CONVERT DAY TO SEC
      CALL os('X=CX    ',x=k2,c=sectoday)
!
!     COMPUTE RS (DONE IN DIFSOU)
!
!----------------------------------------------------------------------
!     LET'S NOW COMPUTE SOURCE TERMS
!
!     FIRST TRACER O2 (RANK IND_O2)
!     warning: here there are lots of choices that can be changed
!              by the user:
!              1- reareration is considered for only water surface,
!              2- bentic demand concerns only the bottom or all the water
!              column (here the second choice is considered)
!              these choices and many others, have to be decided by WAQ
!              experts
!
!     EXPLICIT PART
!     =============
!
!     SURFACE SOURCES
      DO i=1,npoin2
        j=(nplan-1)*npoin2+i
        texp%ADR(ind_o2)%P%R(j) = texp%ADR(ind_o2)%P%R(j) +
     &    t31%R(j)*k2%R(i)*max((t21%R(i)-tn%ADR(ind_o2)%P%R(j)),0.d0) -
!                 O2 CONCENTRATION CAN NOT BE GREATER THAN SATURATED
     &    t32%R(i)/max(eps,zprop%R(j))
      ENDDO
!
!     ADD PHOTOSYNTHESIS - RESPIRATION
!
      CALL os('X=X+C   ',x=texp%ADR(ind_o2)%P,c=pmoinr)
!
!     THE REMAINING TERMS
!
      CALL os('X=X+CY  ',x=texp%ADR(ind_o2)%P,y=tn%ADR(ind_ol)%P,
     &         c=-conv_k1 )
      CALL os('X=X+CY  ',x=texp%ADR(ind_o2)%P,y=tn%ADR(ind_nh4)%P,
     &         c=-conv_k44)
!     IMPLICIT PART:
!     ==============
!     SOFAR, ALL TERMS ARE EXPLICIT FOR O2. CHANGE
!     IF THERE ARE DISAPPOINTING RESULTS
!
!     SECOND TRACER [L] ORGANIC LOAD
!
!      CALL OS('X=C     ',X=TIMP%ADR(IND_OL)%P,C=CONV_K1)
      CALL os('X=X+C   ',x=timp%ADR(ind_ol)%P,c=conv_k1)
!
!     THIRD TRACER [NH4]
!
!      CALL OS('X=C     ',X=TIMP%ADR(IND_NH4)%P,C=CONV_K44)
      CALL os('X=X+C   ',x=timp%ADR(ind_nh4)%P,c=conv_k44)
!
!     MASS BALANCE: MASS ADDED BY EXPLICIT TERMS (TO CHECK)
!
!
!-----------------------------------------------------------------------
!
      RETURN
      END