C $Header: /u/gcmpack/MITgcm/pkg/dic/car_flux.F,v 1.16 2010/08/18 12:20:13 mlosch Exp $
C $Name:  $

#include "DIC_OPTIONS.h"

CBOP
C !ROUTINE: CAR_FLUX

C !INTERFACE: ==========================================================
      SUBROUTINE CAR_FLUX( CAR_S, cflux,
     I           bi,bj,imin,imax,jmin,jmax,
     I           myIter,myTime,myThid)

C !DESCRIPTION:
C Calculate carbonate fluxes

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIC_VARS.h"

C !INPUT PARAMETERS: ===================================================
C  myThid               :: thread number
C  myIter               :: current timestep
C  myTime               :: current time
C  CAR_S                :: carbonate source
      INTEGER myIter
      _RL myTime
      INTEGER myThid
      _RL  CAR_S(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER imin, imax, jmin, jmax, bi, bj

C !OUTPUT PARAMETERS: ===================================================
C  cflux                :: carbonate flux
      _RL cflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)

#if (defined ALLOW_PTRACERS && defined DIC_BIOTIC)

C !LOCAL VARIABLES: ====================================================
C  i,j,k                  :: loop indices
c  ko                     :: loop-within-loop index
c  caexport               :: flux of carbonate from base each "productive"
c                            layer
c  depth_l                :: depth and lower interface
c  flux_u, flux_l         :: flux through upper and lower interfaces
c  reminFac               :: abbreviation
c  zbase                  :: depth of bottom of current productive layer
      INTEGER I,J,k, ko, kop1
      _RL zbase
      _RL caexport(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL reminFac
      _RL depth_l
      _RL flux_u  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL flux_l
CEOP

C- Calculate carbonate flux from base of each layer
      DO k=1,nlev
       DO j=jmin,jmax
        DO i=imin,imax
         caexport(i,j) = 0. _d 0
         IF (hFacC(i,j,k,bi,bj).GT.0. _d 0) THEN
C--   If no layer below initial layer (because of bottom or
C--   topography), then remineralize in here
          IF (k.EQ.Nr) THEN
           cflux(i,j,k)=cflux(i,j,k) + CAR_S(i,j,k) 
          ELSEIF ( _hFacC(i,j,k+1,bi,bj).EQ.0. _d 0 ) THEN
           cflux(i,j,k)=cflux(i,j,k) + CAR_S(i,j,k) 
          ELSE
C- flux out of layer k
           caexport(i,j) = CAR_S(i,j,k)*drF(k) * _hFacC(i,j,k,bi,bj) 
          ENDIF
         ENDIF
        ENDDO
       ENDDO
C--   If availabe, flux carbon export downward;
C--   calculate flux to each layer from base of k
       zbase=-rF(k+1)
C--   Upper flux
       DO j=jmin,jmax
        DO i=imin,imax
         flux_u(i,j)  = caexport(i,j)
        ENDDO
       ENDDO
C     Instead of running the loop to ko=Nr and masking the last
C     flux_l, let ko reach only Nr-1 and do a special loop for ko=Nr,
C     in order to save a few expensive exp-function calls
       DO ko=k+1,Nr-1
        kop1   = MIN(Nr,ko+1)
#ifndef NONLIN_FRSURF
C     For the linear free surface, hFacC can be omitted, buying another
C     performance increase of a factor of six on a vector computer.
C     For now this is not implemented via run time flags, in order to 
C     avoid making this code too complicated.
        depth_l  = -rF(ko) + drF(ko)
        reminFac = exp(-(depth_l-zbase)/zca)
#endif /* NONLIN_FRSURF */
        DO j=jmin,jmax
         DO i=imin,imax
          IF ( caexport(i,j) .NE. 0. _d 0 ) THEN
C--   Lower flux (no flux to ocean bottom)
#ifdef NONLIN_FRSURF
           depth_l  = -rF(ko) + drF(ko) * _hFacC(i,j,ko,bi,bj)
           reminFac = exp(-(depth_l-zbase)/zca)
#endif /* NONLIN_FRSURF */
           flux_l   = caexport(i,j)*reminFac
     &          *maskC(i,j,kop1,bi,bj)
C           
           cflux(i,j,ko)=cflux(i,j,ko) + (flux_u(i,j)-flux_l)
     &          *recip_drF(ko) * _recip_hFacC(i,j,ko,bi,bj)
C--   Store flux through upper layer for the next k-level
           flux_u(i,j)  = flux_l
C     endif carexport .ne. 0
          ENDIF
C     i,j-loops
         ENDDO
        ENDDO
C     ko-loop
       ENDDO
C     now do ko = Nr
       ko = Nr
       flux_l = 0. _d 0
       DO j=jmin,jmax
        DO i=imin,imax
         cflux(i,j,ko)=cflux(i,j,ko) + (flux_u(i,j)-flux_l)
     &        *recip_drF(ko) * _recip_hFacC(i,j,ko,bi,bj)
        ENDDO
       ENDDO
C     k-loop
      ENDDO
c
#endif /* defined ALLOW_PTRACERS && defined DIC_BIOTIC */
      RETURN
      END