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

#include "DIC_OPTIONS.h"

CBOP
C !ROUTINE: PHOS_FLUX

C !INTERFACE: ==========================================================
      SUBROUTINE PHOS_FLUX( BIOac, pflux, exportflux,
     I           bi,bj,imin,imax,jmin,jmax,
     I           myIter,myTime,myThid)

C !DESCRIPTION:
C Calculate the PO4 flux to depth from bio activity

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  BIOac                :: biological productivity
      INTEGER myIter
      _RL myTime
      INTEGER myThid
      _RL  BIOac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER imin, imax, jmin, jmax, bi, bj

C !OUTPUT PARAMETERS: ===================================================
C  pflux                :: changes to PO4 due to flux and reminerlization
      _RL pflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL exportflux(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  bexport                :: flux of phosphorus 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 bexport(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 PO4 flux from base of each layer
      DO k=1,nlev
       DO j=jmin,jmax
        DO i=imin,imax
         bexport(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
           pflux(i,j,k)=pflux(i,j,k)+BIOac(i,j,k)*(1. _d 0-DOPfraction) 
          ELSEIF (hFacC(i,j,k+1,bi,bj).EQ.0. _d 0) THEN
           pflux(i,j,k)=pflux(i,j,k)+BIOac(i,j,k)*(1. _d 0-DOPfraction)
          ELSE
C- flux out of layer k
           bexport(i,j)=BIOac(i,j,k)*(1. _d 0-DOPfraction)
     &          *drF(k) * _hFacC(i,j,k,bi,bj)
          ENDIF
         ENDIF
        ENDDO
       ENDDO
C--   If available, flux phosphate 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)  = bexport(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)
C       reminFac = (depth_l/zbase)**(-Kremin)
C     The following form does the same, but is faster 
        reminFac = exp(-Kremin*log(depth_l/zbase))
#endif
        DO j=jmin,jmax
         DO i=imin,imax
          IF ( bexport(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)
C          reminFac = (depth_l/zbase)**(-Kremin)
C     The following form does the same, but is faster 
           reminFac = exp(-Kremin*log(depth_l/zbase))
#endif
           flux_l   = bexport(i,j)*reminFac
     &          *maskC(i,j,kop1,bi,bj)
C     
           pflux(i,j,ko)=pflux(i,j,ko) + (flux_u(i,j)-flux_l)
     &          *recip_drF(ko) * _recip_hFacC(i,j,ko,bi,bj)
           exportflux(i,j,ko)=exportflux(i,j,ko)+flux_u(i,j)
C--   Store flux through upper layer for the next k-level
           flux_u(i,j) = flux_l
C     endif bexport .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
         pflux(i,j,ko)=pflux(i,j,ko) + (flux_u(i,j)-flux_l)
     &        *recip_drF(ko) * _recip_hFacC(i,j,ko,bi,bj)
         exportflux(i,j,ko)=exportflux(i,j,ko)+flux_u(i,j)
        ENDDO
       ENDDO
C     k-loop
      ENDDO
c
#endif /* defined ALLOW_PTRACERS && defined DIC_BIOTIC */
      RETURN
      END