!! Copyright (C) Stichting Deltares, 2012-2014. !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License version 3, !! as published by the Free Software Foundation. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !! contact: delft3d.support@deltares.nl !! Stichting Deltares !! P.O. Box 177 !! 2600 MH Delft, The Netherlands !! !! All indications and logos of, and references to registered trademarks !! of Stichting Deltares remain the property of Stichting Deltares. All !! rights reserved. subroutine methox ( pmsa , fl , ipoint , increm , noseg , & noflux , iexpnt , iknmrk , noq1 , noq2 , & noq3 , noq4 ) !>\file !> Oxidation of methane (New and generic) C********************************************************************** C +----------------------------------------+ C | D E L F T H Y D R A U L I C S | C | WAter Resources and Environment | C +----------------------------------------+ C C*********************************************************************** C C Project : SLIK, ONTW. BODEM-WATER UITWISSELINGSMODULES, Q2935.30 C Author : Johannes Smits C Date : 020524 Version : 0.01 C C History : C C Date Author Description C ------ -------------- ----------------------------------- C 020524 Johannes Smits New generic process for methane oxidation C 090107 Boderie Add light limitation (Guerin) C C*********************************************************************** C C Description of the module : C C General water quality module for DELWAQ: C Methane oxidation kinetics composed of zeroth order and C MM-kinetics for methane and dissolved oxygen or sulphate. C Process is valid for overlying water as well as sediment. C C ----- description of parameters ----- C Name T L I/O Description Units C ---- --- - - ------------------- ---- C CCH4 R*4 1 I concentration of dissolved methane [gC/m3] C CHFUNC R*4 1 L function for CH4 effect on the oxidation rate [-] C COX R*4 1 I concentration of dissolved oxygen [g/m3] C COXC R*4 1 I critical oxygen conc. for methane oxidation [g/m3] C CRTEMP R*4 1 I critical temperature for nitrification [oC] C CSU R*4 1 I concentration of sulphate [gS/m3] C CSUC R*4 1 I critical sulphate conc. for methane oxidation [gS/m3] C DELT R*4 1 I timestep for processes [d] C FL (1) R*4 1 O methane oxidation flux with DO [gC/m3/d] C FL (2) R*4 1 O methane oxidation flux with sulphate [gC/m3/d] C K0OXI1 R*4 1 I zeroth order methane oxidation rate for DO [gC/m3/d] C K0OXI2 R*4 1 I zeroth order methane oxid. rate for sulphate [gC/m3/d] C KOXI1 R*4 1 I MM methane oxidation rate for DO [gC/m3/d] C KOXI2 R*4 1 I MM methane oxidation rate for sulphate [gC/m3/d] C KSCH4 R*4 1 I half saturation constant for methane [gC/m3] C KSOX R*4 1 I half saturation constant for oxygen [g/m3] C KSSU R*4 1 I half saturation constant for sulphate [gS/m3] C KTOXI1 R*4 1 I temperature coefficient for oxidation with DO [-] C KTOXI2 R*4 1 I temperature coefficient for oxidation with sulphate [-] C KSRadFr R*4 1 I half saturation for light inhibition for fraction of C surface radiance (-) C KSRadSh R*4 1 I shape factor for light inhibition function (-) C OXFUNC R*4 1 O function for DO effect on the oxidation rate [-] C POROS R*4 1 L porosity [-] C RadInh R*4 1 I irradiation level with 100% inhibition (W/m2) C Rad R*4 1 I irradiation at the segment upper-boundary (W/m2) C SUFUNC R*4 1 O function for sulphate effect on the oxidation rate [-] C TEMP R*4 1 I ambient temperature [oC] C TEMPC R*4 1 L ambient temperature correction function [-] C TEMP20 R*4 1 L ambient temperature - stand. temp (20) [oC] C C Logical Units : - C Modules called : - C Name Type Library C ------ ----- ------------ C IMPLICIT NONE C REAL PMSA ( * ) , FL (*) INTEGER IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX, + IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4 C INTEGER IP1, IP2, IP3, IP4, IP5, IP6, IP7, IP8, IP9, IP10, + IP11, IP12, IP13, IP14, IP15, IP16, IP17, IP18, IP19, IP20, + IP21, IP22, IP23, IP24, IP25 INTEGER IN1, IN2, IN3, IN4, IN5, IN6, IN7, IN8, IN9, IN10, + IN11, IN12, IN13, IN14, IN15, IN16, IN17, IN18, IN19, IN20, + IN21, IN22, IN23, IN24, IN25 INTEGER IFLUX , ISEG , IKMRK1, ILUMON C REAL CCH4 , COX , CSU REAL K0OXI1 , K0OXI2 , KOXI1 , KOXI2 , KSCH4 , KSOX , + KSSU , COXC , CSUC , CHFUNC , OXFUNC , SUFUNC , + LIFUNC REAL POROS , CRTEMP , KTOXI1 , KTOXI2 , TEMP , TEMPC , + TEMP20 , KSRADFR, KSRADSH REAL DELT, RAD, RADINH, RADFR50INH REAL COX_RATIO REAL CSU_RATIO REAL FLCOX REAL FLCSU LOGICAL FIRST SAVE FIRST DATA FIRST /.TRUE./ C CALL GETMLU(ILUMON) C C IN1 = INCREM( 1) IN2 = INCREM( 2) IN3 = INCREM( 3) IN4 = INCREM( 4) IN5 = INCREM( 5) IN6 = INCREM( 6) IN7 = INCREM( 7) IN8 = INCREM( 8) IN9 = INCREM( 9) IN10 = INCREM(10) IN11 = INCREM(11) IN12 = INCREM(12) IN13 = INCREM(13) IN14 = INCREM(14) IN15 = INCREM(15) IN16 = INCREM(16) IN17 = INCREM(17) IN18 = INCREM(18) IN19 = INCREM(19) IN20 = INCREM(20) IN21 = INCREM(21) IN22 = INCREM(22) IN23 = INCREM(23) IN24 = INCREM(24) IN25 = INCREM(25) C IP1 = IPOINT( 1) IP2 = IPOINT( 2) IP3 = IPOINT( 3) IP4 = IPOINT( 4) IP5 = IPOINT( 5) IP6 = IPOINT( 6) IP7 = IPOINT( 7) IP8 = IPOINT( 8) IP9 = IPOINT( 9) IP10 = IPOINT(10) IP11 = IPOINT(11) IP12 = IPOINT(12) IP13 = IPOINT(13) IP14 = IPOINT(14) IP15 = IPOINT(15) IP16 = IPOINT(16) IP17 = IPOINT(17) IP18 = IPOINT(18) IP19 = IPOINT(19) IP20 = IPOINT(20) IP21 = IPOINT(21) IP22 = IPOINT(22) IP23 = IPOINT(23) IP24 = IPOINT(24) IP25 = IPOINT(25) C C -----Warnings----- C IF (FIRST) THEN IF (PMSA(IP7).LE. 0.0) THEN WRITE (ILUMON, *) 'WARNING : half saturation constant', + ' KsMet should be greater than zero' ELSEIF (PMSA(IP8).LE. 0.0) THEN WRITE (ILUMON, *) 'WARNING : half saturation constant', + ' KsOxMet should be greater than zero' ELSEIF (PMSA(IP13) .LE. 0.0) THEN WRITE (ILUMON, *) 'WARNING : half saturation constant', + ' KsSuMet should be greater than zero' ELSEIF (PMSA(IP17) .LE. 0.0) THEN WRITE (ILUMON, *) 'WARNING : Poros should be greater', + 'than zero' ENDIF FIRST = .FALSE. ENDIF C IFLUX = 0 DO 9000 ISEG = 1 , NOSEG !! CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1) !! IF ( IKMRK1 .GT. 0) THEN IF (BTEST(IKNMRK(ISEG),0)) THEN C CCH4 = MAX ( 0.0, PMSA(IP1 ) ) COX = MAX ( 0.0, PMSA(IP2 ) ) CSU = MAX ( 0.0, PMSA(IP3 ) ) K0OXI1 = PMSA(IP4 ) KOXI1 = PMSA(IP5 ) KTOXI1 = PMSA(IP6 ) KSCH4 = PMSA(IP7 ) KSOX = PMSA(IP8 ) COXC = PMSA(IP9 ) K0OXI2 = PMSA(IP10) KOXI2 = PMSA(IP11) KTOXI2 = PMSA(IP12) KSSU = PMSA(IP13) CSUC = PMSA(IP14) TEMP = PMSA(IP15) CRTEMP = PMSA(IP16) POROS = PMSA(IP17) DELT = PMSA(IP18) RADINH= PMSA(IP19) RAD = PMSA(IP20) KSRADFr= PMSA(IP21) KSRADSh= PMSA(IP22) C C Set the rates according to CRTEMP, COXC and CSUC C IF ( TEMP .LT. CRTEMP .OR. + COX .LE. 0.0 .OR. COX .LE. (COXC * POROS) ) THEN KOXI1 = 0.0 ENDIF IF ( COX .LE. 0.0 .OR. COX .GT. (COXC * POROS) ) THEN K0OXI1 = 0.0 ENDIF C IF ( TEMP .LT. CRTEMP .OR. + CSU .LE. 0.0 .OR. CSU .LE. (CSUC * POROS) ) THEN KOXI2 = 0.0 ENDIF IF ( CSU .LE. 0.0 .OR. CSU .GT. (CSUC * POROS) ) THEN K0OXI2 = 0.0 ENDIF C IF ( COX .GT. (COXC * POROS) ) THEN KOXI2 = 0.0 K0OXI2 = 0.0 ENDIF C C Calculate both methane oxidation fluxes C TEMP20 = TEMP - 20.0 TEMPC = KTOXI1 ** TEMP20 C IF ( (KSCH4 * POROS + CCH4) .GT. 0.0) THEN CHFUNC = CCH4 / ( KSCH4 * POROS + CCH4 ) ELSE CHFUNC = 0.0 ENDIF C IF ( (KSOX * POROS + COX) .GT. 0.0) THEN OXFUNC = COX / ( KSOX * POROS + COX ) ELSE OXFUNC = 0.0 ENDIF C IF ( (KSSU * POROS + CSU) .GT. 0.0) THEN SUFUNC = CSU / ( KSSU * POROS + CSU ) ELSE SUFUNC = 0.0 ENDIF C C Light inhibition function - analogue to salinity dependant mort function IF ( KSRADFR .GE. 0.0 .AND. KSRADFR .LE. 1.0 + .AND. RADINH .GT. 0.1) THEN LIFUNC =-1/(1 + EXP(KSRADSh*(min(RAD/RADINH,1.0) - KSRADFR)))+ 1 ELSE LIFUNC = 0.0 ENDIF C FLCOX = K0OXI1 + KOXI1 * TEMPC * CHFUNC * OXFUNC * (1-LIFUNC) FLCSU = K0OXI2 + KOXI2 * TEMPC * CHFUNC * SUFUNC * (1-LIFUNC) COX_RATIO = 5.33 FLCOX = MIN(FLCOX,0.5*COX/COX_RATIO/DELT) FLCOX = MIN(FLCOX,0.9*CCH4/DELT) CSU_RATIO = 2.67 FLCSU = MIN(FLCSU,0.9*CSU/CSU_RATIO/DELT) FLCSU = MIN(FLCSU,0.9*CCH4/DELT) FL( 1+IFLUX ) = FLCOX FL( 2+IFLUX ) = FLCSU C C Oxygen and sulphate functions are output C PMSA(IP23) = OXFUNC PMSA(IP24) = SUFUNC PMSA(IP25) = LIFUNC C ENDIF C IFLUX = IFLUX + NOFLUX IP1 = IP1 + IN1 IP2 = IP2 + IN2 IP3 = IP3 + IN3 IP4 = IP4 + IN4 IP5 = IP5 + IN5 IP6 = IP6 + IN6 IP7 = IP7 + IN7 IP8 = IP8 + IN8 IP9 = IP9 + IN9 IP10 = IP10 + IN10 IP11 = IP11 + IN11 IP12 = IP12 + IN12 IP13 = IP13 + IN13 IP14 = IP14 + IN14 IP15 = IP15 + IN15 IP16 = IP16 + IN16 IP17 = IP17 + IN17 IP18 = IP18 + IN18 IP19 = IP19 + IN19 IP20 = IP20 + IN20 IP21 = IP21 + IN21 IP22 = IP22 + IN22 IP23 = IP23 + IN23 IP24 = IP24 + IN24 IP25 = IP25 + IN25 C 9000 CONTINUE C RETURN C END