calcs2d_o2.f

Go to the documentation of this file.

!                     *********************
                      SUBROUTINE calcs2d_o2
!                     *********************
!
     & (npoin,wattemp,o2satu,demben,formk2,k1,k44,k22,
     &  photo,resp,tn,texp,timp,t2,t3,t4,hprop,un,vn,debug)
!
!***********************************************************************
! WAQTEL   V8P1
!***********************************************************************
!
!brief    COMPUTES SOURCE TERMS FOR O2 WAQ PROCESS
!
!history  R. ATA
!+        21/09/2014
!+        V7P0
!+       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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DEBUG         |-->| FOR DEBUGGING
!| DEMBEN        |-->| BENTIC DEMAND
!| FORMK2        |-->| FORMULA FOR COMPUTINK K2
!| K1            |-->| CONST. OF DEGRADATION KINETIC OF ORGANIC LOAD
!| HPROP         |-->| PROPAGATION DEPTH
!| K22           |<--| CONST. OF REAERATION
!| K44           |-->| CONST. OF NITRIFICATION KINETIC
!| NPOIN         |-->| NUMBER OF NODES IN THE MESH
!| O2SATU        |<--| O2 SATURATION DENSITY OF WATER (CS)
!| PHOTO         |-->| PHOTOSYNTHESIS
!| WATTEMP       |-->| TEMPERATURE OF WATER (DEG CELSIUS)
!| T2,T3,T4      |-->| WORKING ARRAY
!| TEXP          |<--| EXPLICIT SOURCE TERM.
!| TIMP          |<--| IMPLICIT SOURCE TERM.
!| TN            |-->| TRACERS AT TIME N
!| UN,VN         |-->| VELOCITY COMPONENTS AT TIME N
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

      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_calcs2d_o2 => calcs2d_o2
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      INTEGER          , INTENT(IN   ) :: formk2,npoin,debug
      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
      TYPE(bief_obj)   , INTENT(INOUT) :: texp,timp,t2,t3,t4
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     LOCAL VARIABLES
      INTEGER                     :: i
      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
!
!-----------------------------------------------------------------------
!
!
!     PRELIMINARY COMPUTATIONS
!
      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,npoin
        IF( ind_t.GT.0 ) power=tn%ADR(ind_t)%P%R(i)-20.d0
!       CORR2T AND BENCOR STOCKED HERE IN T3,T4
        t3%R(i)=corr2**power
        t4%R(i)=bencorr*(corr1**power)
      ENDDO
!
!     COMPUTE CS (O2SATU, STOCKED IN T2)
!
      IF( ind_t.EQ.0 )THEN
        CALL satur_o2(o2satu,formcs,wattemp,eps)
        CALL os('X=C     ',x=t2,c=o2satu       )
      ELSE
        DO i=1,npoin
          CALL satur_o2(t2%R(i),formcs,tn%ADR(ind_t)%P%R(i),eps)
        ENDDO
      ENDIF
!
!     COMPUTE K2
!
      CALL reaer(formk2,k2,k22,npoin,1,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,npoin
        texp%ADR(ind_o2)%P%R(i) = texp%ADR(ind_o2)%P%R(i)
     &    + t3%R(i)*k2%R(i) * max((t2%R(i)-tn%ADR(ind_o2)%P%R(i)),0.d0)
!                    O2 CONCENTRATION CAN NOT BE GREATER THAN SATURATED
     &    - t4%R(i)/max(eps,hprop%R(i))
      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)
      IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP1:',texp%ADR(ind_o2)%P%R(1)

!     IMPLICIT PART:
!     ==============
!     SOFAR, ALL TERMS ARE EXPLICIT FOR O2. CHANGE
!     IF THERE ARE DISAPPOINTING RESULTS
!
!     SECOND TRACER [L] ORGANIC LOAD
!
! /!\ contrary to 3D formulation
      CALL os('X=X+CY  ',x=timp%ADR(ind_ol)%P,y=hprop,c=-conv_k1)
!      CALL OS('X=X+C   ',X=TIMP%ADR(IND_OL)%P,C=CONV_K1)
      IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP2:',timp%ADR(ind_ol)%P%R(1)
!
!     THIRD TRACER [NH4]
!
! /!\ contrary to 3D formulation
      CALL os('X=X+CY  ',x=timp%ADR(ind_nh4)%P,y=hprop,c=-conv_k44)
!      CALL OS('X=X+C   ',X=TIMP%ADR(IND_NH4)%P,C=CONV_K44)
      IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP3:',timp%ADR(ind_nh4)%P%R(1)
!
!-----------------------------------------------------------------------
!
      RETURN
      END