New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 775 for branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90 – NEMO

Ignore:
Timestamp:
2007-12-19T14:45:15+01:00 (16 years ago)
Author:
gm
Message:

dev_001_GM - PISCES in F90 : encapsulation of all p4z...F files in module F90 + doctor norme for local variables - compilation OK

File:
1 moved

Legend:

Unmodified
Added
Removed
  • branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90

    r774 r775  
    1  
    2 CCC $Header$  
    3 CCC  TOP 1.0 , LOCEAN-IPSL (2005)  
    4 C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    5 C --------------------------------------------------------------------------- 
    6 CDIR$ LIST 
    7       SUBROUTINE p4zche 
    8 #if defined key_top && defined key_pisces 
    9 CCC--------------------------------------------------------------------- 
    10 CCC 
    11 CCC          ROUTINE p4zche : PISCES MODEL 
    12 CCC          ***************************** 
    13 CCC 
    14 CCC     PURPOSE. 
    15 CCC     -------- 
    16 CCC          *P4ZCHE* : Sea water chemistry computed following OCMIP protocol 
    17 CCC 
    18 CCC 
    19 CC     EXTERNALS. 
    20 CC     ---------- 
    21 CC          rhop 
    22 CC 
    23 CC   MODIFICATIONS: 
    24 CC   -------------- 
    25 CC      original :      1988 E. Maier-Reimer 
    26 CC      additions :     1998 O. Aumont 
    27 CC      modifications : 1999 C. Le Quere 
    28 CC      modifications : 2004 O. Aumont 
    29 CC      modifications : 2006 R. Gangsto 
    30 CC---------------------------------------------------------------------- 
    31 CC parameters and commons 
    32 CC ====================== 
    33 CDIR$ nolist 
    34       USE oce_trc 
    35       USE trp_trc 
    36       USE sms 
    37       IMPLICIT NONE 
     1MODULE p4zche 
     2   !!====================================================================== 
     3   !!                         ***  MODULE p4zche  *** 
     4   !! TOP :   PISCES Sea water chemistry computed following OCMIP protocol 
     5   !!====================================================================== 
     6   !! History :    -   !  1988     (E. Maier-Reimer)  Original code 
     7   !!              -   !  1998     (O. Aumont)  addition 
     8   !!              -   !  1999     (C. Le Quere)  modification 
     9   !!             1.0  !  2004     (O. Aumont)  modification 
     10   !!              -   !  2006     (R. Gangsto)  modification 
     11   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
     12   !!---------------------------------------------------------------------- 
     13#if defined key_pisces 
     14   !!---------------------------------------------------------------------- 
     15   !!   'key_pisces'                                       PISCES bio-model 
     16   !!---------------------------------------------------------------------- 
     17   !!   p4z_che        :  Sea water chemistry computed following OCMIP protocol 
     18   !!---------------------------------------------------------------------- 
     19   USE oce_trc         ! 
     20   USE trp_trc         ! 
     21   USE sms             !  
     22 
     23   IMPLICIT NONE 
     24   PRIVATE 
     25 
     26   PUBLIC   p4z_che    ! called in p4zprg.F90 
     27 
     28   !!* Substitution 
    3829#include "domzgr_substitute.h90" 
    39 CDIR$ list 
    40 CC---------------------------------------------------------------------- 
    41 CC local declarations 
    42 CC ================== 
    43 C 
    44       INTEGER ji, jj, jk 
    45       REAL tkel, sal,  qtt, zbuf1, zbuf2 
    46       REAL pres, tc, cl, cpexp, cek0, oxy, cpexp2 
    47       REAL zsqrt, ztr, zlogt, cek1 
    48       REAL zqtt, qtt2, sal15, zis, zis2, zisqrt 
    49       REAL ckb, ck1, ck2, ckw, ak1, ak2, akb, aksp0, akw 
    50       REAL ckp1, ckp2, ckp3, cksi, akp1, akp2, akp3, aksi 
    51       REAL st, ft, cks, ckf, aks, akf, aksp1 
    52  
    53 C 
    54 C* 1. CHEMICAL CONSTANTS - SURFACE LAYER 
    55 C --------------------------------------- 
    56 C 
    57       DO jj = 1,jpj 
    58         DO ji = 1,jpi 
    59 C 
    60 C* 1.1 SET ABSOLUTE TEMPERATURE 
    61 C ------------------------------ 
    62 C 
    63           tkel = tn(ji,jj,1)+273.16 
    64           qtt = tkel*0.01 
    65           qtt2=qtt*qtt 
    66           sal = sn(ji,jj,1) + (1.-tmask(ji,jj,1))*35. 
    67           zqtt=log(qtt) 
    68 C 
    69 C* 1.2 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    70 C      AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    71 C ------------------------------------------------------- 
    72 C 
    73           cek0 = c00+c01/qtt+c02*zqtt+sal*(c03+c04*qtt+c05*qtt2) 
    74           cek1 = ca0+ca1/qtt+ca2*zqtt+ca3*qtt2+sal*(ca4 
    75      &      +ca5*qtt+ca6*qtt2) 
    76 C 
    77 C* 1.3 LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970) 
    78 C ------------------------------------------------------------ 
    79 C 
    80           oxy = ox0+ox1/qtt+ox2*zqtt+sal*(ox3+ox4*qtt+ox5*qtt2) 
    81 C 
    82 C* 1.4 SET SOLUBILITIES OF O2 AND CO2 
    83 C ----------------------------------- 
    84 C 
    85           chemc(ji,jj,1) = exp(cek0)*1.E-6*rhop(ji,jj,1)/1000. 
    86           chemc(ji,jj,2) = exp(oxy)*oxyco 
    87           chemc(ji,jj,3) = exp(cek1)*1.E-6*rhop(ji,jj,1)/1000. 
    88 C 
    89         ENDDO 
     30   !!---------------------------------------------------------------------- 
     31   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)  
     32   !! $Header:$  
     33   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     34   !!---------------------------------------------------------------------- 
     35 
     36CONTAINS 
     37 
     38   SUBROUTINE p4z_che 
     39      !!--------------------------------------------------------------------- 
     40      !!                     ***  ROUTINE p4z_che  *** 
     41      !! 
     42      !! ** Purpose :   Sea water chemistry computed following OCMIP protocol 
     43      !! 
     44      !! ** Method  : - ... 
     45      !!--------------------------------------------------------------------- 
     46      INTEGER  ::   ji, jj, jk 
     47      REAL(wp) ::   ztkel, zsal , zqtt  , zbuf1 , zbuf2 
     48      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zcek0, zoxy  , zcpexp2 
     49      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1 
     50      REAL(wp) ::   zlqtt, zqtt2, zsal15, zis   , zis2 , zisqrt 
     51      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
     52      REAL(wp) ::   zckp1, zckp2, zckp3 , zcksi , zakp1, zakp2 , zakp3, zaksi 
     53      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaks , zakf  , zaksp1 
     54      !!--------------------------------------------------------------------- 
     55 
     56      ! CHEMICAL CONSTANTS - SURFACE LAYER 
     57      ! ---------------------------------- 
     58 
     59      DO jj = 1, jpj 
     60         DO ji = 1, jpi 
     61 
     62            !                             ! SET ABSOLUTE TEMPERATURE 
     63            ztkel = tn(ji,jj,1) + 273.16 
     64            zqtt  = ztkel * 0.01 
     65            zqtt2 = zqtt * zqtt 
     66            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35. 
     67            zlqtt = LOG( zqtt ) 
     68 
     69            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
     70            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
     71            zcek0 = c00 + c01 / zqtt + c02 * zlqtt + zsal * ( c03 + c04 * zqtt + c05 * zqtt2 ) 
     72            zcek1 = ca0 + ca1 / zqtt + ca2 * zlqtt + ca3 * zqtt2 + zsal*( ca4 + ca5 * zqtt + ca6 * zqtt2 ) 
     73 
     74            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970) 
     75            zoxy  = ox0 + ox1 / zqtt + ox2 * zlqtt + zsal * ( ox3 + ox4 * zqtt + ox5 * zqtt2 ) 
     76 
     77            !                             ! SET SOLUBILITIES OF O2 AND CO2 
     78            chemc(ji,jj,1) = EXP( zcek0 ) * 1.e-6 * rhop(ji,jj,1) / 1000. 
     79            chemc(ji,jj,2) = EXP( zoxy  ) * oxyco 
     80            chemc(ji,jj,3) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. 
     81 
     82         END DO 
    9083      END DO 
    91 C 
    92 C* 2 CHEMICAL CONSTANTS - DEEP OCEAN 
    93 C ------------------------------------- 
    94 C 
    95       DO jk = 1,jpk 
    96         DO jj = 1,jpj 
    97           DO ji = 1,jpi 
    98 C 
    99 C* 2.1 SET PRESSION 
    100 C ----------------- 
    101 C 
    102             pres = 1.025e-1*fsdept(ji,jj,jk) 
    103 C 
    104 C* 2.2 SET ABSOLUTE TEMPERATURE 
    105 C ------------------------------ 
    106 C 
    107             tkel   = tn(ji,jj,jk)+273.16 
    108             qtt    = tkel*0.01 
    109             sal    = sn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*35. 
    110             zsqrt  = sqrt(sal) 
    111             sal15  = zsqrt*sal 
    112             zlogt  = log(tkel) 
    113             ztr    = 1./tkel 
    114             zis    = 19.924*sal/(1000.-1.005*sal) 
    115             zis2   = zis*zis 
    116             zisqrt = sqrt(zis) 
    117             tc = tn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*20. 
    118 C 
    119 C* 2.3 CHLORINITY (WOOSTER ET AL., 1969) 
    120 C --------------------------------------- 
    121 C 
    122             cl = sal*salchl 
    123 C 
    124 C* 2.4 TOTAL SULFATE CONCENTR. [MOLES/kg soln] 
    125 C -------------------------------------------- 
    126 C 
    127             st = st1*cl*st2 
    128 C 
    129 C* 2.5 TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 
    130 C --------------------------------------------- 
    131 C 
    132             ft = ft1*cl*ft2 
    133 C 
    134 C* 2.6 DISSOCIATION CONSTANT FOR SULFATES 
    135 C on free H scale (Dickson 1990) 
    136 C ------------------------------------------------------- 
    137 C 
    138             cks=exp(ks1*ztr+ks0+ks2*zlogt+(ks3*ztr+ks4+ks5*zlogt) 
    139      &      *zisqrt+(ks6*ztr+ks7+ks8*zlogt)*zis+ks9*ztr*zis 
    140      &      *zisqrt+ks10*ztr*zis2+log(ks11+ks12*sal)) 
    141 C 
    142 C* 2.7 DISSOCIATION CONSTANT FOR FLUORIDES 
    143 C on free H scale (Dickson and Riley 79) 
    144 C ------------------------------------------------------- 
    145 C 
    146             ckf=exp(kf1*ztr+kf0+kf2*zisqrt+log(kf3+kf4*sal)) 
    147  
    148 C 
    149 C* 2.4 DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
    150 C ------------------------------------------------------- 
    151 C 
    152             ckb = (cb0+cb1*zsqrt+cb2*sal+cb3*sal15+cb4*sal*sal)*ztr 
    153      &          +(cb5+cb6*zsqrt+cb7*sal)+ 
    154      &          (cb8+cb9*zsqrt+cb10*sal)*zlogt+cb11*zsqrt*tkel 
    155      &          +log((1.+st/cks+ft/ckf)/(1.+st/cks)) 
    156             ck1 = c10*ztr+c11+c12*zlogt+c13*sal+c14*sal**2 
    157             ck2 = c20*ztr+c21+c22*sal+c23*sal**2 
    158 C 
    159 C* 2.5 PKW (H2O) (DICKSON AND RILEY, 1979) 
    160 C ----------------------------------------- 
    161 C 
    162             ckw = cw0*ztr+cw1+cw2*zlogt+(cw3*ztr+cw4+cw5*zlogt)* 
    163      &          zsqrt+cw6*sal 
    164  
    165 C 
    166 C 
    167 C* 2.10 DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale) 
    168 C  --------------------------------------------------------------------- 
    169 C 
    170             ckp1 = cp10+cp11*ztr+cp12*zlogt+zsqrt*(cp13*ztr 
    171      &           +cp14)+sal*(cp15*ztr+cp16) 
    172             ckp2 = cp20+cp21*ztr+cp22*zlogt+zsqrt*(cp23*ztr 
    173      &           +cp24)+sal*(cp25*ztr+cp26) 
    174             ckp3 = cp30+cp31*ztr+zsqrt*(cp32*ztr 
    175      &           +cp33)+sal*(cp34*ztr+cp35) 
    176             cksi = cs10+cs11*ztr+cs12*zlogt+zisqrt*(cs13*ztr 
    177      &           +cs14)+zis*(cs15*ztr+cs16)+zis2*(cs17*ztr 
    178      &           +cs18)+log(1.+cs19*sal) 
    179      &           +log(cs20+cs21*sal) 
    180  
    181 C 
    182 C*2.7 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
    183 C       (S=27-43, T=2-25 DEG C) AT pres =0 (ATMOSPH. PRESSURE) 
    184 C       (MUCCI 1983) 
    185 C ------------------------------------------------------------- 
    186 C 
    187             aksp0 = akcc1+akcc2*tkel+akcc3*ztr+akcc4*log10(tkel)+ 
    188      &      (akcc5+akcc6*tkel+ 
    189      &      akcc7*ztr)*zsqrt+akcc8*sal+akcc9*sal15 
    190  
    191  
    192 C 
    193 C* 2.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
    194 C ----------------------------------------------------------------- 
    195 C 
    196             ak1   = 10**(ck1) 
    197             ak2   = 10**(ck2) 
    198             akb   = exp(ckb) 
    199             akp1  = exp(ckp1) 
    200             akp2  = exp(ckp2) 
    201             akp3  = exp(ckp3) 
    202             aksi  = exp(cksi) 
    203             akw   = exp(ckw) 
    204             aksp1 = 10**(aksp0) 
    205             aks   = exp(cks) 
    206             akf   = exp(ckf) 
    207  
    208  
    209 C 
    210 C* 2.8 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970) 
    211 C        (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE 
    212 C        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 
    213 C        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN 
    214 C        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 
    215 C        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 
    216 C        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 
    217 C        AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL 
    218 C        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 
    219 C        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON 
    220 C        P. 1285)) 
    221 C ----------------------------------------------------------- 
    222 C 
    223             cpexp = pres /(rgas*tkel) 
    224             cpexp2 = pres * pres/(rgas*tkel) 
    225 C 
    226 C* 2.9 KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 
    227 C        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 
    228 C        (CF. BROECKER ET AL., 1982) 
    229 C -------------------------------------------------------- 
    230 C 
    231             zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc) 
    232             zbuf2 = 0.5*(devk4(3)+devk5(3)*tc) 
    233             akb3(ji,jj,jk) = akb*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    234  
    235             zbuf1 = -(devk1(1)+devk2(1)*tc+devk3(1)*tc*tc) 
    236             zbuf2 = 0.5*(devk4(1)+devk5(1)*tc) 
    237             ak13(ji,jj,jk) = ak1*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    238  
    239             zbuf1 = -(devk1(2)+devk2(2)*tc+devk3(2)*tc*tc) 
    240             zbuf2 = 0.5*(devk4(2)+devk5(2)*tc) 
    241             ak23(ji,jj,jk) = ak2*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    242  
    243             zbuf1 = -(devk1(4)+devk2(4)*tc+devk3(4)*tc*tc) 
    244             zbuf2 = 0.5*(devk4(4)+devk5(4)*tc) 
    245             akp13(ji,jj,jk) = akp1*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    246  
    247             zbuf1 = -(devk1(5)+devk2(5)*tc+devk3(5)*tc*tc) 
    248             zbuf2 = 0.5*(devk4(5)+devk5(5)*tc) 
    249             akp23(ji,jj,jk) = akp2*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    250  
    251             zbuf1 = -(devk1(6)+devk2(6)*tc+devk3(6)*tc*tc) 
    252             zbuf2 = 0.5*(devk4(6)+devk5(6)*tc) 
    253             akp33(ji,jj,jk) = akp3*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    254  
    255             zbuf1 = -(devk1(7)+devk2(7)*tc+devk3(7)*tc*tc) 
    256             zbuf2 = 0.5*(devk4(7)+devk5(7)*tc) 
    257             akw3(ji,jj,jk) = akw*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    258  
    259 C  Ksi 
    260 C            aksi3(ji,jj,jk) = aksi 
    261 C 
    262 C  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys) 
    263 C  "deltaVsi and deltaKsi have been estimated from the value of boric acid" 
    264 C 
    265             zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc) 
    266             zbuf2 = 0.5*(devk4(3)+devk5(3)*tc) 
    267             aksi3(ji,jj,jk) = aksi*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    268  
    269 C 
    270 C 
    271 C* 2.15 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE  
    272 C        AS FUNCTION OF PRESSURE FOLLOWING MILLERO 
    273 C        (P. 1285) AND BERNER (1976) 
    274 C ------------------------------------------------- 
    275  
    276             zbuf1 = -(devk1(8)+devk2(8)*tc+devk3(8)*tc*tc) 
    277             zbuf2 = 0.5*(devk4(8)+devk5(8)*tc) 
    278             aksp(ji,jj,jk) = aksp1*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    279  
    280 C  Pressure correction for sulfate and fluoride 
    281 C 
    282             zbuf1 = -(devk1(9)+devk2(9)*tc+devk3(9)*tc*tc) 
    283             zbuf2 = 0.5*(devk4(9)+devk5(9)*tc) 
    284             aks3(ji,jj,jk) = aks*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    285  
    286             zbuf1 = -(devk1(10)+devk2(10)*tc+devk3(10)*tc*tc) 
    287             zbuf2 = 0.5*(devk4(10)+devk5(10)*tc) 
    288             akf3(ji,jj,jk) = akf*exp(zbuf1*cpexp+zbuf2*cpexp2) 
    289  
    290  
    291 C 
    292 C* 2.11 TOTAL BORATE CONCENTR. [MOLES/L] 
    293 C -------------------------------------- 
    294 C 
    295             borat(ji,jj,jk) = bor1*cl*bor2 
    296 C 
    297 C  2.12 Iron and SIO3 saturation concentration from ... 
    298 C  ---------------------------------------------------- 
    299 C 
    300          sio3eq(ji,jj,jk)=exp(log(10.)*(6.44-968./tkel))*1E-6 
    301          fekeq(ji,jj,jk)=10**(17.27-1565.7/(273.15+tc)) 
    302 C 
    303           ENDDO 
    304         ENDDO 
     84 
     85      ! CHEMICAL CONSTANTS - DEEP OCEAN 
     86      ! ------------------------------- 
     87 
     88      DO jk = 1, jpk 
     89         DO jj = 1, jpj 
     90            DO ji = 1, jpi 
     91 
     92               ! SET PRESSION 
     93               zpres   = 1.025e-1 * fsdept(ji,jj,jk) 
     94 
     95               ! SET ABSOLUTE TEMPERATURE 
     96               ztkel   = tn(ji,jj,jk) + 273.16 
     97               zqtt    = ztkel * 0.01 
     98               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 
     99               zsqrt  = SQRT( zsal ) 
     100               zsal15  = zsqrt * zsal 
     101               zlogt  = LOG( ztkel ) 
     102               ztr    = 1. / ztkel 
     103               zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal ) 
     104               zis2   = zis * zis 
     105               zisqrt = SQRT( zis ) 
     106               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 
     107 
     108               ! CHLORINITY (WOOSTER ET AL., 1969) 
     109               zcl     = zsal * salchl 
     110 
     111               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 
     112               zst     = st1 * zcl * st2 
     113 
     114               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 
     115               zft     = ft1 * zcl * ft2 
     116 
     117               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
     118               zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           & 
     119                  &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   & 
     120                  &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      & 
     121                  &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  ) 
     122 
     123               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
     124               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  ) 
     125 
     126               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
     127               zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   & 
     128                  &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            & 
     129                  &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             & 
     130                  &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  ) 
     131!!gm zsal**2 to be replaced by a *... 
     132               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal**2 
     133               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2 
     134 
     135               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
     136               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal 
     137 
     138               ! DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale) 
     139               zckp1 = cp10 + cp11 * ztr + cp12 * zlogt + zsqrt * ( cp13 * ztr + cp14 ) + zsal * ( cp15 * ztr + cp16 ) 
     140               zckp2 = cp20 + cp21 * ztr + cp22 * zlogt + zsqrt * ( cp23 * ztr + cp24 ) + zsal * ( cp25 * ztr + cp26 ) 
     141               zckp3 = cp30 + cp31 * ztr                + zsqrt * ( cp32 * ztr + cp33 ) + zsal * ( cp34 * ztr + cp35 ) 
     142               zcksi = cs10 + cs11 * ztr + cs12 * zlogt + zisqrt* ( cs13 * ztr + cs14 ) + zis * ( cs15 * ztr + cs16 )  & 
     143                  &                                    + zis2  * ( cs17 * ztr + cs18 ) + LOG( 1.   + cs19 * zsal )     & 
     144                  &                                                                    + LOG( cs20 + cs21 * zsal ) 
     145 
     146               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
     147               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 
     148               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   & 
     149                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15 
     150 
     151               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
     152               zak1    = 10**(zck1) 
     153               zak2    = 10**(zck2) 
     154               zakb    = EXP( zckb  ) 
     155               zakp1   = EXP( zckp1 ) 
     156               zakp2   = EXP( zckp2 ) 
     157               zakp3   = EXP( zckp3 ) 
     158               zaksi   = EXP( zcksi ) 
     159               zakw    = EXP( zckw ) 
     160               zaksp1  = 10**(zaksp0) 
     161               zaks    = exp( zcks ) 
     162               zakf    = exp( zckf ) 
     163 
     164               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 
     165               !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE 
     166               !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 
     167               !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN 
     168               !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 
     169               !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 
     170               !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 
     171               !        & GIESKES (1970), P. 1285-1286 (THE SMALL 
     172               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 
     173               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 
     174               zcpexp  = zpres /(rgas*ztkel) 
     175               zcpexp2 = zpres * zpres/(rgas*ztkel) 
     176 
     177               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 
     178               !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 
     179               !        (CF. BROECKER ET AL., 1982) 
     180               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc ) 
     181               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc ) 
     182               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     183 
     184               zbuf1  = -(devk1(1)+devk2(1)*ztc+devk3(1)*ztc*ztc) 
     185               zbuf2  = 0.5*(devk4(1)+devk5(1)*ztc) 
     186               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     187 
     188               zbuf1  =     - ( devk1(2) + devk2(2) * ztc + devk3(2) * ztc * ztc ) 
     189               zbuf2  = 0.5 * ( devk4(2) + devk5(2) * ztc ) 
     190               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     191 
     192               zbuf1  =     - ( devk1(4) + devk2(4) * ztc + devk3(4) * ztc * ztc ) 
     193               zbuf2  = 0.5 * ( devk4(4) + devk5(4) * ztc ) 
     194               akp13(ji,jj,jk) = zakp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     195 
     196               zbuf1  =     - ( devk1(5) + devk2(5) * ztc + devk3(5) * ztc * ztc ) 
     197               zbuf2  = 0.5 * ( devk4(5) + devk5(5) * ztc ) 
     198               akp23(ji,jj,jk) = zakp2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     199 
     200               zbuf1  =     - ( devk1(6) + devk2(6) * ztc + devk3(6) * ztc * ztc ) 
     201               zbuf2  = 0.5 * ( devk4(6) + devk5(6) * ztc ) 
     202               akp33(ji,jj,jk) = zakp3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     203 
     204               zbuf1  =     - ( devk1(7) + devk2(7) * ztc + devk3(7) * ztc * ztc ) 
     205               zbuf2  = 0.5 * ( devk4(7) + devk5(7) * ztc ) 
     206               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     207 
     208               !  Ksi 
     209               !            aksi3(ji,jj,jk) = zaksi 
     210               ! 
     211               !  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys) 
     212               !  "deltaVsi and deltaKsi have been estimated from the value of boric acid" 
     213               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc ) 
     214               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc ) 
     215               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     216 
     217               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE  
     218               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO 
     219               !        (P. 1285) AND BERNER (1976) 
     220               zbuf1  =     - ( devk1(8) + devk2(8) * ztc + devk3(8) * ztc * ztc ) 
     221               zbuf2  = 0.5 * ( devk4(8) + devk5(8) * ztc ) 
     222               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     223 
     224               !  Pressure correction for sulfate and fluoride 
     225               zbuf1  =     - ( devk1(9) + devk2(9) * ztc + devk3(9) * ztc * ztc ) 
     226               zbuf2  = 0.5 * ( devk4(9) + devk5(9) * ztc ) 
     227               aks3(ji,jj,jk) = zaks   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     228 
     229               zbuf1  =     - ( devk1(10) + devk2(10) * ztc + devk3(10) * ztc * ztc ) 
     230               zbuf2  = 0.5 * ( devk4(10) + devk5(10) * ztc ) 
     231               akf3(ji,jj,jk) = zakf   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     232 
     233               ! TOTAL BORATE CONCENTR. [MOLES/L] 
     234               borat(ji,jj,jk) = bor1 * zcl * bor2 
     235 
     236               ! Iron and SIO3 saturation concentration from ... 
     237               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6 
     238               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) ) 
     239 
     240            END DO 
     241         END DO 
    305242      END DO 
    306 C 
    307 #endif 
    308 C 
    309       RETURN 
    310       END 
     243      ! 
     244   END SUBROUTINE p4z_che 
     245 
     246#else 
     247   !!====================================================================== 
     248   !!  Dummy module :                                   No PISCES bio-model 
     249   !!====================================================================== 
     250CONTAINS 
     251   SUBROUTINE p4z_che( kt )                   ! Empty routine 
     252      INTEGER, INTENT( in ) ::   kt 
     253      WRITE(*,*) 'p4z_che: You should not have seen this print! error?', kt 
     254   END SUBROUTINE p4z_che 
     255#endif  
     256 
     257   !!====================================================================== 
     258END MODULE  p4zche 
Note: See TracChangeset for help on using the changeset viewer.