#include "BLING_OPTIONS.h"

CBOP
      subroutine BLING_CARBONATE_SYS(
     I           PTR_DIC, PTR_ALK, PTR_PO4,
#ifdef USE_SIBLING
     I           PTR_SI,
#endif
     I           bi, bj, imin, imax, jmin, jmax,
     I           myTime, myIter, myThid)

C     =================================================================
C     | subroutine bling_carbonate_sys
C     | o Calculate carbonate fluxes
C     |   Also update pH (3d field)
C     =================================================================

      IMPLICIT NONE

C     == GLobal variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "BLING_VARS.h"

C     == Routine arguments ==
C     PTR_DIC              :: dissolved inorganic carbon
C     PTR_ALK              :: alkalinity
C     PTR_PO4              :: phosphate
C     myThid               :: thread Id. number
C     myIter               :: current timestep
C     myTime               :: current time
      INTEGER myThid
      INTEGER myIter
      _RL  myTime
      _RL  PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL  PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL  PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef USE_SIBLING
      _RL  PTR_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#endif
      INTEGER imin, imax, jmin, jmax, bi, bj

#ifdef ALLOW_PTRACERS

C     == Local variables ==
C     i,j,k             :: loop indices
C     carbonate         :: local value of calcium carbonate
C     calcium           :: local value of Ca
C     diclocal          :: local value of DIC
C     alklocal          :: local value of ALK
C     pCO2local         :: local value of pCO2
C     pHlocal           :: local value of pH
C     CO3iter           :: iterations counter for CO3 ion calculation
C     CO3iterMax        :: total number of iterations
C     silicaDEEP        :: subsurface silica concentration
       INTEGER i,j,k
       _RL carbonate
       _RL calcium
       _RL po4local
       _RL diclocal
       _RL alklocal
       _RL pCO2local
       _RL pHlocal
       _RL silicaDEEP

       _RL ttmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL stmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

c      INTEGER CO3iter
c      INTEGER CO3iterMax
CEOP

C  Since pH is now a 3D field and is solved for at every time step
C  few iterations are needed
c      CO3iterMax = 1

C determine carbonate ion concentration through full domain
C determine calcite saturation state

C$TAF LOOP = parallel
       DO k=1,Nr

          DO j=jMin,jMax
            DO i=iMin,iMax
             ttmp(i,j) = theta(i,j,k,bi,bj)
             stmp(i,j) = salt(i,j,k,bi,bj)
            ENDDO
          ENDDO

#ifdef CARBONCHEM_SOLVESAPHE
cav  I'm using ttmp and stmp instead of surftemp and surfsalt like
cav  Jonathan / they are not surf when k>1
cav  Does new carbon_chem in DICpkg need theta (3d) or ttmp (2d)???

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_CALL(
     &   'DIC_COEFFS_DEEP',myThid)
#endif
C Calculate carbon coefficients
        CALL DIC_COEFFS_SURF(
     I                       ttmp,stmp,
     I                       bi,bj,iMin,iMax,jMin,jMax,myThid)

C Now correct the coefficients for pressure dependence
        CALL DIC_COEFFS_DEEP(
     I                       ttmp,stmp,
     I                       bi,bj,iMin,iMax,jMin,jMax,
     I                       k,myThid)

#else

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_CALL(
     &   'CARBON_COEFFS_PRESSURE_DEP',myThid)
#endif

C  Get coefficients for carbonate calculations
          CALL CARBON_COEFFS_PRESSURE_DEP(
     I                       ttmp, stmp,
     I                       bi, bj, imin, imax, jmin, jmax,
     I                       k, myThid)
#endif

C--------------------------------------------------

C$TAF LOOP = parallel
          DO j=jMin,jMax
C$TAF LOOP = parallel
           DO i=iMin,iMax

             IF ( hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
C$TAF init dic_caco3 = static, 2

#ifdef CARBONCHEM_SOLVESAPHE
             calcium = cat(i,j,bi,bj)
#else
C  Estimate calcium concentration from salinity
             calcium = 1.028 _d -2*salt(i,j,k,bi,bj)/35. _d 0
#endif

             po4local = PTR_PO4(i,j,k)
             diclocal = PTR_DIC(i,j,k)
             alklocal = PTR_ALK(i,j,k)
             pHlocal  = pH(i,j,k,bi,bj)

C  Assume constant deep silica value
C  30 micromol = 0.03 mol m-3
C  unless SIBLING is used
#ifdef USE_SIBLING
             silicaDEEP = PTR_SI(i,j,k)
#else
             silicaDEEP = 0.03 _d 0
#endif

#ifdef CARBONCHEM_SOLVESAPHE
      IF ( selectPHsolver.GT.0 ) THEN
C Use Munhoven (2013) Solvesaphe routine to calculate pH and pCO2
#ifdef ALLOW_DEBUG
         IF (debugMode) CALL DEBUG_CALL('AHINI_FOR_AT',myThid)
#endif
CAV since we carry pH, no need for an initial guess
C call AHINI_FOR_AT to get better initial guess of pH
c               CALL AHINI_FOR_AT(alklocal*permil,
c     I           diclocal*permil,
c     I           bt(i,j,bi,bj),
c     O           pHlocal,
c     I           i,j,k,bi,bj,myIter,myThid )
#ifdef ALLOW_DEBUG
         IF (debugMode) CALL DEBUG_CALL('CALC_PCO2_SOLVESAPHE',myThid)
#endif
            CALL CALC_PCO2_SOLVESAPHE(
     I          theta(i,j,k,bi,bj),salt(i,j,k,bi,bj),
     I          diclocal, po4local,
     I          silicaDEEP,alklocal,
     U          pHlocal,pCO2local,carbonate,
     I          i,j,k,bi,bj,myIter,myThid )
      ELSE
C Use the original Follows et al. (2006) solver
#endif /* CARBONCHEM_SOLVESAPHE */
#ifdef ALLOW_DEBUG
            IF (debugMode) CALL DEBUG_CALL('CALC_PCO2_APPROX',myThid)
#endif

C  Evaluate carbonate (CO3) ions concentration
C  iteratively

c             DO CO3iter = 1, CO3iterMax

C--------------------------------------------------

               CALL CALC_PCO2_APPROX(
     I          ttmp(i,j),stmp(i,j),
     I          diclocal, po4local,
     I          silicaDEEP,alklocal,
     I          ak1(i,j,bi,bj),ak2(i,j,bi,bj),
     I          ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
     I          aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
     I          aksi(i,j,bi,bj),akf(i,j,bi,bj),
     I          ak0(i,j,bi,bj),fugf(i,j,bi,bj),ff(i,j,bi,bj),
     I          bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
     U          pHlocal,pCO2local,carbonate,
     I          i,j,k,bi,bj,myIter,myThid )
c             ENDDO

#ifdef CARBONCHEM_SOLVESAPHE
      ENDIF
#endif /* CARBONCHEM_SOLVESAPHE */

              pH(i,j,k,bi,bj) = pHlocal

C  Calculate calcium carbonate (calcite and aragonite)
C  saturation state
             omegaC(i,j,k,bi,bj) = calcium * carbonate /
     &                          Ksp_TP_Calc(i,j,bi,bj)
             omegaAr(i,j,k,bi,bj) = calcium * carbonate /
     &                          Ksp_TP_Arag(i,j,bi,bj)

           else

             pH(i,j,k,bi,bj) = 0. _d 0
             omegaC(i,j,k,bi,bj)  = 0. _d 0
             omegaAr(i,j,k,bi,bj) = 0. _d 0

           endif

         ENDDO
        ENDDO
       ENDDO

#endif /* ALLOW_PTRACERS */
       RETURN
       END
