Changeset 617 for trunk/NEMO/TOP_SRC/SMS/p4zche.F

Timestamp:

2007-02-21T14:03:10+01:00 (17 years ago)

Author:

opalod

Message:

* empty log message *

File:

• trunk/NEMO/TOP_SRC/SMS/p4zche.F (modified) (8 diffs)

trunk/NEMO/TOP_SRC/SMS/p4zche.F

@@ -14 +14 @@
 CCC     PURPOSE.
 CCC     --------
-CCC          *P4ZCHE* SETS CHEMICAL CONSTANTS 
+CCC          *P4ZCHE* : Sea water chemistry computed following OCMIP protocol
 CCC
 CCC
@@ -27 +27 @@
 CC      modifications : 1999 C. Le Quere
 CC      modifications : 2004 O. Aumont
+CC      modifications : 2006 R. Gangsto
 CC----------------------------------------------------------------------
 CC parameters and commons
@@ -45 +46 @@
       REAL pres, tc, cl, cpexp, cek0, oxy, cpexp2
       REAL zsqrt, ztr, zlogt, cek1
-      REAL zqtt, qtt2, sal15, zis, zis2
+      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
@@ -110 +114 @@
             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
@@ -116 +121 @@
 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
@@ -123 +153 @@
      &          +(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
@@ -131 +162 @@
             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.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
 C -----------------------------------------------------------------
 C
-            ak1 = 10**(ck1)
-            ak2 = 10**(ck2)
-            akb = exp(ckb)
-            akw = exp(ckw)
+            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.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       (INGLE, 1800, EQ. 6)
+C       (MUCCI 1983)
 C -------------------------------------------------------------
 C
-            aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)
-     &          +akcc4*tkel*tkel)
+            aksp0 = akcc1+akcc2*tkel+akcc3*ztr+akcc4*log10(tkel)+
+     &      (akcc5+akcc6*tkel+
+     &      akcc7*ztr)*zsqrt+akcc8*sal+akcc9*sal15
+
 C
 C* 2.8 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
@@ -183 +242 @@
             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
-C  2.10 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE (OR ARAGONITE)
-C        AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970)
+
+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 -----------------------------------------------------------------
-C
-            aksp(ji,jj,jk) = aksp0*exp(cpexp*(devks-devkst*tc))
+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]