CCC $Header$ 
CCC  TOP 1.0 , LOCEAN-IPSL (2005) 
C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
C ---------------------------------------------------------------------------
CDIR$ LIST
      SUBROUTINE p4zche
#if defined key_passivetrc && defined key_trc_pisces
CCC---------------------------------------------------------------------
CCC
CCC          ROUTINE p4zche : PISCES MODEL
CCC          *****************************
CCC
CCC     PURPOSE.
CCC     --------
CCC          *P4ZCHE* : Sea water chemistry computed following OCMIP protocol
CCC
CCC
CC     EXTERNALS.
CC     ----------
CC          rhop
CC
CC   MODIFICATIONS:
CC   --------------
CC      original :      1988 E. Maier-Reimer
CC      additions :     1998 O. Aumont
CC      modifications : 1999 C. Le Quere
CC      modifications : 2004 O. Aumont
CC      modifications : 2006 R. Gangsto
CC----------------------------------------------------------------------
CC parameters and commons
CC ======================
CDIR$ nolist
      USE oce_trc
      USE trp_trc
      USE sms
      IMPLICIT NONE
#include "domzgr_substitute.h90"
CDIR$ list
CC----------------------------------------------------------------------
CC local declarations
CC ==================
C
      INTEGER ji, jj, jk
      REAL tkel, sal,  qtt, zbuf1, zbuf2
      REAL pres, tc, cl, cpexp, cek0, oxy, cpexp2
      REAL zsqrt, ztr, zlogt, cek1
      REAL zqtt, qtt2, sal15, zis, zis2, zisqrt
      REAL ckb, ck1, ck2, ckw, ak1, ak2, akb, aksp0, akw
      REAL ckp1, ckp2, ckp3, cksi, akp1, akp2, akp3, aksi
      REAL st, ft, cks, ckf, aks, akf, aksp1

C
C* 1. CHEMICAL CONSTANTS - SURFACE LAYER
C ---------------------------------------
C
      DO jj = 1,jpj
        DO ji = 1,jpi
C
C* 1.1 SET ABSOLUTE TEMPERATURE
C ------------------------------
C
          tkel = tn(ji,jj,1)+273.16
          qtt = tkel*0.01
          qtt2=qtt*qtt
          sal = sn(ji,jj,1) + (1.-tmask(ji,jj,1))*35.
          zqtt=log(qtt)
C
C* 1.2 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
C      AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
C -------------------------------------------------------
C
          cek0 = c00+c01/qtt+c02*zqtt+sal*(c03+c04*qtt+c05*qtt2)
          cek1 = ca0+ca1/qtt+ca2*zqtt+ca3*qtt2+sal*(ca4
     &      +ca5*qtt+ca6*qtt2)
C
C* 1.3 LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970)
C ------------------------------------------------------------
C
          oxy = ox0+ox1/qtt+ox2*zqtt+sal*(ox3+ox4*qtt+ox5*qtt2)
C
C* 1.4 SET SOLUBILITIES OF O2 AND CO2
C -----------------------------------
C
          chemc(ji,jj,1) = exp(cek0)*1.E-6*rhop(ji,jj,1)/1000.
          chemc(ji,jj,2) = exp(oxy)*oxyco
          chemc(ji,jj,3) = exp(cek1)*1.E-6*rhop(ji,jj,1)/1000.
C
        ENDDO
      END DO
C
C* 2 CHEMICAL CONSTANTS - DEEP OCEAN
C -------------------------------------
C
      DO jk = 1,jpk
        DO jj = 1,jpj
          DO ji = 1,jpi
C
C* 2.1 SET PRESSION
C -----------------
C
            pres = 1.025e-1*fsdept(ji,jj,jk)
C
C* 2.2 SET ABSOLUTE TEMPERATURE
C ------------------------------
C
            tkel   = tn(ji,jj,jk)+273.16
            qtt    = tkel*0.01
            sal    = sn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*35.
            zsqrt  = sqrt(sal)
            sal15  = zsqrt*sal
            zlogt  = log(tkel)
            ztr    = 1./tkel
            zis    = 19.924*sal/(1000.-1.005*sal)
            zis2   = zis*zis
            zisqrt = sqrt(zis)
            tc = tn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*20.
C
C* 2.3 CHLORINITY (WOOSTER ET AL., 1969)
C ---------------------------------------
C
            cl = sal*salchl
C
C* 2.4 TOTAL SULFATE CONCENTR. [MOLES/kg soln]
C --------------------------------------------
C
            st = st1*cl*st2
C
C* 2.5 TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
C ---------------------------------------------
C
            ft = ft1*cl*ft2
C
C* 2.6 DISSOCIATION CONSTANT FOR SULFATES
C on free H scale (Dickson 1990)
C -------------------------------------------------------
C
            cks=exp(ks1*ztr+ks0+ks2*zlogt+(ks3*ztr+ks4+ks5*zlogt)
     &      *zisqrt+(ks6*ztr+ks7+ks8*zlogt)*zis+ks9*ztr*zis
     &      *zisqrt+ks10*ztr*zis2+log(ks11+ks12*sal))
C
C* 2.7 DISSOCIATION CONSTANT FOR FLUORIDES
C on free H scale (Dickson and Riley 79)
C -------------------------------------------------------
C
            ckf=exp(kf1*ztr+kf0+kf2*zisqrt+log(kf3+kf4*sal))

C
C* 2.4 DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
C -------------------------------------------------------
C
            ckb = (cb0+cb1*zsqrt+cb2*sal+cb3*sal15+cb4*sal*sal)*ztr
     &          +(cb5+cb6*zsqrt+cb7*sal)+
     &          (cb8+cb9*zsqrt+cb10*sal)*zlogt+cb11*zsqrt*tkel
     &          +log((1.+st/cks+ft/ckf)/(1.+st/cks))
            ck1 = c10*ztr+c11+c12*zlogt+c13*sal+c14*sal**2
            ck2 = c20*ztr+c21+c22*sal+c23*sal**2
C
C* 2.5 PKW (H2O) (DICKSON AND RILEY, 1979)
C -----------------------------------------
C
            ckw = cw0*ztr+cw1+cw2*zlogt+(cw3*ztr+cw4+cw5*zlogt)*
     &          zsqrt+cw6*sal

C
C
C* 2.10 DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale)
C  ---------------------------------------------------------------------
C
            ckp1 = cp10+cp11*ztr+cp12*zlogt+zsqrt*(cp13*ztr
     &           +cp14)+sal*(cp15*ztr+cp16)
            ckp2 = cp20+cp21*ztr+cp22*zlogt+zsqrt*(cp23*ztr
     &           +cp24)+sal*(cp25*ztr+cp26)
            ckp3 = cp30+cp31*ztr+zsqrt*(cp32*ztr
     &           +cp33)+sal*(cp34*ztr+cp35)
            cksi = cs10+cs11*ztr+cs12*zlogt+zisqrt*(cs13*ztr
     &           +cs14)+zis*(cs15*ztr+cs16)+zis2*(cs17*ztr
     &           +cs18)+log(1.+cs19*sal)
     &           +log(cs20+cs21*sal)

C
C*2.7 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
C       (S=27-43, T=2-25 DEG C) AT pres =0 (ATMOSPH. PRESSURE)
C       (MUCCI 1983)
C -------------------------------------------------------------
C
            aksp0 = akcc1+akcc2*tkel+akcc3*ztr+akcc4*log10(tkel)+
     &      (akcc5+akcc6*tkel+
     &      akcc7*ztr)*zsqrt+akcc8*sal+akcc9*sal15


C
C* 2.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
C -----------------------------------------------------------------
C
            ak1   = 10**(ck1)
            ak2   = 10**(ck2)
            akb   = exp(ckb)
            akp1  = exp(ckp1)
            akp2  = exp(ckp2)
            akp3  = exp(ckp3)
            aksi  = exp(cksi)
            akw   = exp(ckw)
            aksp1 = 10**(aksp0)
            aks   = exp(cks)
            akf   = exp(ckf)


C
C* 2.8 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
C        (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE
C        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
C        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
C        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
C        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
C        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
C        AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL
C        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
C        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON
C        P. 1285))
C -----------------------------------------------------------
C
            cpexp = pres /(rgas*tkel)
            cpexp2 = pres * pres/(rgas*tkel)
C
C* 2.9 KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
C        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
C        (CF. BROECKER ET AL., 1982)
C --------------------------------------------------------
C
            zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc)
            zbuf2 = 0.5*(devk4(3)+devk5(3)*tc)
            akb3(ji,jj,jk) = akb*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(1)+devk2(1)*tc+devk3(1)*tc*tc)
            zbuf2 = 0.5*(devk4(1)+devk5(1)*tc)
            ak13(ji,jj,jk) = ak1*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(2)+devk2(2)*tc+devk3(2)*tc*tc)
            zbuf2 = 0.5*(devk4(2)+devk5(2)*tc)
            ak23(ji,jj,jk) = ak2*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(4)+devk2(4)*tc+devk3(4)*tc*tc)
            zbuf2 = 0.5*(devk4(4)+devk5(4)*tc)
            akp13(ji,jj,jk) = akp1*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(5)+devk2(5)*tc+devk3(5)*tc*tc)
            zbuf2 = 0.5*(devk4(5)+devk5(5)*tc)
            akp23(ji,jj,jk) = akp2*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(6)+devk2(6)*tc+devk3(6)*tc*tc)
            zbuf2 = 0.5*(devk4(6)+devk5(6)*tc)
            akp33(ji,jj,jk) = akp3*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(7)+devk2(7)*tc+devk3(7)*tc*tc)
            zbuf2 = 0.5*(devk4(7)+devk5(7)*tc)
            akw3(ji,jj,jk) = akw*exp(zbuf1*cpexp+zbuf2*cpexp2)

C  Ksi
C            aksi3(ji,jj,jk) = aksi
C
C  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys)
C  "deltaVsi and deltaKsi have been estimated from the value of boric acid"
C
            zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc)
            zbuf2 = 0.5*(devk4(3)+devk5(3)*tc)
            aksi3(ji,jj,jk) = aksi*exp(zbuf1*cpexp+zbuf2*cpexp2)

C
C
C* 2.15 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
C        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
C        (P. 1285) AND BERNER (1976)
C -------------------------------------------------

            zbuf1 = -(devk1(8)+devk2(8)*tc+devk3(8)*tc*tc)
            zbuf2 = 0.5*(devk4(8)+devk5(8)*tc)
            aksp(ji,jj,jk) = aksp1*exp(zbuf1*cpexp+zbuf2*cpexp2)

C  Pressure correction for sulfate and fluoride
C
            zbuf1 = -(devk1(9)+devk2(9)*tc+devk3(9)*tc*tc)
            zbuf2 = 0.5*(devk4(9)+devk5(9)*tc)
            aks3(ji,jj,jk) = aks*exp(zbuf1*cpexp+zbuf2*cpexp2)

            zbuf1 = -(devk1(10)+devk2(10)*tc+devk3(10)*tc*tc)
            zbuf2 = 0.5*(devk4(10)+devk5(10)*tc)
            akf3(ji,jj,jk) = akf*exp(zbuf1*cpexp+zbuf2*cpexp2)


C
C* 2.11 TOTAL BORATE CONCENTR. [MOLES/L]
C --------------------------------------
C
            borat(ji,jj,jk) = bor1*cl*bor2
C
C  2.12 Iron and SIO3 saturation concentration from ...
C  ----------------------------------------------------
C
         sio3eq(ji,jj,jk)=exp(log(10.)*(6.44-968./tkel))*1E-6
         fekeq(ji,jj,jk)=10**(17.27-1565.7/(273.15+tc))
C
          ENDDO
        ENDDO
      END DO
C
#endif
C
      RETURN
      END