C######################################################################
C Turbulence-closure physics, pelagic biology and 
C diagenetic model of the ocean margin (OMEX project MAS2-CT93-0069)
C By:
C Karline Soetaert, Peter M.J. Herman and Jack J Middelburg
C NIOO-CEMO - Vierstraat 28-4401 EA Yerseke, the Netherlands
C Claire Burren, Paul Tett
C UCNW, Bangor Menai Bridge, Gwynedd LL59 5EY, U.K.
C
C File: CARSED.FOR
C This file contains the subroutines CarSed and CarComp
C that describe the two carbon fractions in the sediment
C
C Include files are:
C XCOMMON.DCS: the dimensions of the arrays
C XCOMMON.BLK: the parameter and variable values and names as declared in
C             'INPUT.PAR' and 'INPUT.VAR' 
C XCOMMON.ALL: constants/arrays specific to the model

C######################################################################  

      SUBROUTINE CARSED()
*--------------------------------------------------------------------*
* Calls the bottom carbon submodel CARComp with the appropriate      *
* conversion factor (TimeFac), that converts the time unit of the    *
* submodel (i.c. days) to the time unit of the entire model (hours). *
*--------------------------------------------------------------------*

*-------------
* Declarations
      IMPLICIT NONE

      INCLUDE 'XCOMMON.DCS'
      INCLUDE 'XCOMMON.BLK'
      INCLUDE 'XCOMMON.ALL'

      DOUBLE PRECISION Day2Hour 
      PARAMETER (Day2Hour = 1.d0/24.d0)

*-----------
* Statements


! The interface between pelagic and benthic module (C, N deposition)
      CALL INTERFACE(AmmConvFlux)

! If the bottom is dynamically modeled, calculate C dynamics, else do nothing
      IF (INT(Diagen) .EQ. 1) THEN
        CALL CARCOMP (Day2Hour, 1)
      ELSE IF (INT (DIAGEN) .EQ. 2) THEN

! SEMI-DYNAMIC, VERTICALLY INTEGRATED BOUNDARY

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Decay constants at the temperature
! Q10 factor for mineralisation rate (Q10Deg); 
      Q10fac = EXP((Temp(1) - 20.)*LOG(Q10deg)/10.)

      rFDET = cFDET20*Q10fac
      rSDET = cSDET20*Q10fac
	Kdissi = CKdissi20*Q10fac

! sediment fast detritus decreases by mineralization and 
! increases by flux
      fDegrad = rFDET*BFDET(1)
      dBFDET(1) = dBFDET(1) +(- fDegrad + FDETflux)*Day2Hour
       
! sediment slow detritus decreases by mineralization and 
! increases by flux
      sDegrad = rSDET*BSDET(1)
      dBSDET(1) = dBSDET(1) +(- sDegrad + SDETflux)*Day2Hour

	SIPROD = Kdissi * BDETSI(1)

	DBDETSI(1) = DBDETSI(1) +(- SIPROD + DETSIFLUX)*Day2Hour
	
! Total nitrogen in the model 
! (pelagic units in mmol/m2; benthic units in nmol/cm2)
      TotalN = TotalNPel + (BFDET(1)*NCrFD + BSDET(1)*NCrSD)*.01
	      
      ENDIF                                      ! DIAGEN
 	
      END                                        ! of the main routine

************************************************************
*       DYNAMICS OF DETRITUS IN THE BOTTOM:
************************************************************

        SUBROUTINE CARComp(TimeFac, Switch)
*-------------------------------------------------------------*
* Calculates the rate of change of the two detrital fractions *
* in the bottom                                               *
* as described in Soetaert et al. (1996)                      *
*  Geochim. Cosmochim. Acta: 1019-1040.                       *
*-------------------------------------------------------------*

*-------------
* Declarations

      IMPLICIT NONE

      INCLUDE 'XCOMMON.DCS'
      INCLUDE 'XCOMMON.BLK'
      INCLUDE 'XCOMMON.ALL'

      DOUBLE PRECISION    
     &            C0, AdsCoeff, One, TimeFac, Flux
      INTEGER     I, Mixtype, Portype, Upbound, adstype, Switch
      EXTERNAL    Transsub

*-----------
* Statements
      One = 1.d0

! The number of sediment slices (XBEN) is in XCOMMON.DCS (XBEN)

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
! TRANSPORT of Fast Decay Carbon 

! Variable dispersion coefficient (bioturbation)
      Mixtype = 2

! Solid substance
      Portype = 2

! Upper boundary condition = flux condition
      Upbound = 2

! Flux
      Flux = FDETflux

! No adsorption
      Adstype = 0
      Adscoeff = 0.

      CALL Transsub(XBEN, BFDET, dBFDET, SedRate,
     &       Mixtype, Biotur, One,
     &       portype, porsurf, pordepth, 
     &       Upbound, C0, Flux,
     &       Adstype, Adscoeff, TimeFac)

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
! TRANSPORT of Slow Decay Carbon 

      Flux = SDETflux

      CALL Transsub(XBEN, BSDET, dBSDET, SedRate,
     &       Mixtype, Biotur, One,
     &       portype, porsurf, pordepth, 
     &       Upbound, C0, Flux,
     &       Adstype, Adscoeff, TimeFac)

! Total flux of carbon, expressed in the model's time unit
      Carflux = (Flux + FDETflux)*TimeFac

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
! DEGRADATION of Carbon

! C is consumed at a first-order rate
! First calculate the degradation rate, based on the temperature 
! assume temperature in the bottom = temperature of the lowest pelagic compartment

! Q10 factor for mineralisation rate (Q10Deg); 
      Q10fac = EXP((Temp(1) - 20.)*LOG(Q10deg)/10.)

! Decay constants at the temperature
      rFDET = cFDET20*Q10fac
      rSDET = cSDET20*Q10fac

      DO I=1, XBEN
        dBFDET(I) = dBFDET(I) - rFDET * BFDET(I) * TimeFac
        dBSDET(I) = dBSDET(I) - rSDET * BSDET(I) * TimeFac
      ENDDO

! Average the mineralization rates in each sediment layer; expressed per day      
! Only if a dynamic run has started (Switch = 1);
! During the steady-state calculation of the carbon profiles, no averaging
! (Switch = -1); Last call to carcomp: switch = 0
      IF (Switch .EQ. 1) THEN
        DO I = 1,XBEN
         FDETMin(I) = FDETMin(I) + 
     &                    rFDET * BFDET(I) * CarbonFac       
         SDETMin(I) = SDETMin(I) +
     &                   rSDET * BSDET(I) * CarbonFac       
        ENDDO
      ELSE IF (Switch .EQ. 0) THEN
        DO I = 1, XBEN
         FDETMin(I) = rFDET * BFDET(I)
         SDETMin(I) = rSDET * BSDET(I)
        ENDDO
      ENDIF

      END                                        ! of CARCOMP
