source: trunk/NEMO/TOP_SRC/SMS/p4zche.F @ 247

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* SETS CHEMICAL CONSTANTS 
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----------------------------------------------------------------------
CC parameters and commons
CC ======================
CDIR$ NOLIST
      USE oce_trc
      USE trp_trc
      USE sms
      IMPLICIT NONE
CDIR$ list
CC----------------------------------------------------------------------
CC local declarations
CC ==================
C
      INTEGER ji, jj, jk
      REAL tkel, sal, rrr, qtt
      REAL pres, tc, cl, cpexp
      REAL akb, temzer, cek0, oxy
      REAL zsqrt, ztr, zlogt
      REAL zqtt, qtt2, sal15
      REAL ckb, ck1, ck2, ckw, ak1, ak2, aksp0
CC----------------------------------------------------------------------
CC statement functions
CC ===================
CDIR$ NOLIST
#include "domzgr_substitute.h90"
CDIR$ LIST
C
C* 1. CHEMICAL CONSTANTS - SURFACE LAYER
C ---------------------------------------
      temzer = 273.16
C
C vertical slab
C =============
C
      DO jj = 1,jpj
        DO ji = 1,jpi
C
C* 1.1 SET ABSOLUTE TEMPERATURE
C ------------------------------
C
          tkel = tn(ji,jj,1)+temzer
          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, 1974)
C -------------------------------------------------------
C
          cek0 = c00+c01/qtt+c02*zqtt+sal*(c03+c04*qtt+c05*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 CHEMICAL CHEMICAL CONSTANTS
C --------------------------------------
C
          chemc(ji,jj,1) = exp(cek0)*1.E-6
C
C* 1.5 O2 SOLUBILITY IN SEAWATER (WEISS, 1970, CF. EQ. 4)
C ---------------------------------------------------------
C
          chemc(ji,jj,3) = exp(oxy)*oxyco
C
        ENDDO
      END DO
C
C* 2 CHEMICAL CONSTANTS - DEEP OCEAN
C -------------------------------------
C
      DO jk = 1,jpk
C
C* 2.1 APPROX. SEAWATER PRESSURE AT U-POINT DEPTH (BAR)
C ------------------------------------------------------
C
        DO jj=1,jpj
          DO ji = 1,jpi
C
C* 2.2 SET [H+] (FIRST GUESS)
C ----------------------------
C
            pres = 1.025e-1*fsdept(ji,jj,jk)
            hi(ji,jj,jk) = 1.E-7
C
C* 2.3 SET ABSOLUTE TEMPERATURE
C ------------------------------
C
            tkel   = tn(ji,jj,jk)+temzer
            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
C
C* 2.4 CHLORINITY (WOOSTER ET AL., 1969)
C ---------------------------------------
C
            cl = sal*salchl
C
C* 2.5 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1974)
C -------------------------------------------------------
C
            cek0 = c00+c01/qtt+c02*log(qtt)+
     &          sal*(c03+c04*qtt+c05*qtt*qtt)
C
C  COEFFICIENT OCMIP 
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
            ck1 = c10*ztr+c11+c12*zlogt+(c13*ztr+c14)*zsqrt+
     $          c15*sal+c16*sal15+log(1.+c17*sal)
            ck2 = c20*ztr+c21+c22*zlogt+(c23*ztr+c24)*zsqrt+c25*sal
     $          +c26*sal15+log(1.+c27*sal)
C
C* 2.6 PKW (H2O) (DICKSON AND RILEY, 1979)
C -----------------------------------------
C
            ckw = cw0*ztr+cw1+cw2*zlogt+(cw3*ztr+cw4+cw5*zlogt)*
     $          zsqrt+cw6*sal
C
C* 2.7 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
C -----------------------------------------------------------------
C
            ak1 = exp(ck1)
            ak2 = exp(ck2)
            akb = exp(ckb)
            akw3(ji,jj,jk) = exp(ckw)
C
C*2.8 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
C       (S=27-43, T=2-25 DEG C) AT pres =0 (ATMOSPH. PRESSURE)
C       (INGLE, 1800, EQ. 6)
C -------------------------------------------------------------
C
            aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)
     &          +akcc4*tkel*tkel)
C
C* 2.9 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)
C
C* 2.10 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
            tc = tn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*20.
            akb3(ji,jj,jk) = akb*exp(cpexp*(devkb-devkbt*tc))
            ak13(ji,jj,jk) = ak1*exp(cpexp*(devk1-devk1t*tc))
            ak23(ji,jj,jk) = ak2*exp(cpexp*(devk2-devk2t*tc))
C
C  2.11 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE (OR ARAGONITE)
C        AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970)
C        (P. 1285) AND BERNER (1976)
C -----------------------------------------------------------------
C
            aksp(ji,jj,jk) = aksp0*exp(cpexp*(devks-devkst*tc))
C
C* 2.12 DENSITY OF SEAWATER AND TOTAL BORATE CONCENTR. [MOLES/L]
C ---------------------------------------------------------------
C
            rrr = rhop(ji,jj,jk)/1000.
            borat(ji,jj,jk) = bor1*rrr*cl*bor2
C
C  2.13 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**(16.27-1565.7/(273.15+tn(ji,jj,jk)))
C
          ENDDO
        ENDDO
      END DO
C      
#endif
C
      RETURN
      END