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

Timestamp:

2005-11-14T13:30:28+01:00 (18 years ago)

Author:

opalod

Message:

nemo_v1_update_027 : CE + RB + CT : update of SMS routines

File:

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

trunk/NEMO/TOP_SRC/SMS/p4zche.F

@@ -1 @@
-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
@@ -29 +25 @@
 CC parameters and commons
 CC ======================
-CDIR$ NOLIST
+CDIR$ nolist
       USE oce_trc
       USE trp_trc
       USE sms
       IMPLICIT NONE
+#include "domzgr_substitute.h90"
 CDIR$ list
 CC----------------------------------------------------------------------
@@ -40 +37 @@
 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
+      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
+      REAL ckb, ck1, ck2, ckw, ak1, ak2, akb, aksp0, akw
 C
 C* 1. CHEMICAL CONSTANTS - SURFACE LAYER
 C ---------------------------------------
-      temzer = 273.16
-C
-C vertical slab
-C =============
 C
       DO jj = 1,jpj
@@ -66 +52 @@
 C ------------------------------
 C
-          tkel = tn(ji,jj,1)+temzer
+          tkel = tn(ji,jj,1)+273.16
           qtt = tkel*0.01
           qtt2=qtt*qtt
@@ -72 +58 @@
           zqtt=log(qtt)
 C
-C* 1.2 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1974)
+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)
@@ -82 +71 @@
           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* 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
@@ -99 +85 @@
 C
       DO jk = 1,jpk
-C
-C* 2.1 APPROX. SEAWATER PRESSURE AT U-POINT DEPTH (BAR)
-C ------------------------------------------------------
-C
-        DO jj=1,jpj
+        DO jj = 1,jpj
           DO ji = 1,jpi
 C
-C* 2.2 SET [H+] (FIRST GUESS)
-C ----------------------------
+C* 2.1 SET PRESSION
+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* 2.2 SET ABSOLUTE TEMPERATURE
 C ------------------------------
 C
-            tkel   = tn(ji,jj,jk)+temzer
+            tkel   = tn(ji,jj,jk)+273.16
             qtt    = tkel*0.01
             sal    = sn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*35.
@@ -122 +103 @@
             zlogt  = log(tkel)
             ztr    = 1./tkel
-C
-C* 2.4 CHLORINITY (WOOSTER ET AL., 1969)
+            zis    = 19.924*sal/(1000.-1.005*sal)
+            zis2   = zis*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.5 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1974)
+C* 2.4 DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
 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)
+     &          +(cb5+cb6*zsqrt+cb7*sal)+
+     &          (cb8+cb9*zsqrt+cb10*sal)*zlogt+cb11*zsqrt*tkel
+            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* 2.7 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+     &          zsqrt+cw6*sal
+C
+C* 2.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
 C -----------------------------------------------------------------
 C
-            ak1 = exp(ck1)
-            ak2 = exp(ck2)
+            ak1 = 10**(ck1)
+            ak2 = 10**(ck2)
             akb = exp(ckb)
-            akw3(ji,jj,jk) = exp(ckw)
-C
-C*2.8 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
+            akw = exp(ckw)
+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 -------------------------------------------------------------
 C
-            aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)
+            aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)
      &          +akcc4*tkel*tkel)
 C
-C* 2.9 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
+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
@@ -181 +157 @@
 C
             cpexp = pres /(rgas*tkel)
-C
-C* 2.10 KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+            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
-            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)
+            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)
+            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        (P. 1285) AND BERNER (1976)
@@ -199 +187 @@
             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* 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**(16.27-1565.7/(273.15+tn(ji,jj,jk)))
+         fekeq(ji,jj,jk)=10**(17.27-1565.7/(273.15+tc))
 C
           ENDDO
         ENDDO
       END DO
-C      
+C
 #endif
 C