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 10322 for NEMO/branches/UKMO/dev_r9950_GO8_package/src/TOP/PISCES/SED/sedchem.F90 – NEMO

Ignore:
Timestamp:
2018-11-16T16:06:47+01:00 (5 years ago)
Author:
davestorkey
Message:

UKMO/dev_r9950_GO8_package: Update to be relative to rev 10321 of NEMO4_beta_mirror branch.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/dev_r9950_GO8_package/src/TOP/PISCES/SED/sedchem.F90

    r9950 r10322  
    11MODULE sedchem 
    22 
    3 #if defined key_sed  
    43   !!====================================================================== 
    54   !!                        ***  Module sedchem  *** 
     
    98   USE sed     ! sediment global variable 
    109   USE sedarr 
     10   USE eosbn2, ONLY : neos 
     11   USE lib_mpp         ! distribued memory computing library 
     12 
     13   IMPLICIT NONE 
     14   PRIVATE 
    1115 
    1216   !! * Accessibility 
    13    PUBLIC sed_chem    
     17   PUBLIC sed_chem 
     18   PUBLIC ahini_for_at_sed     ! 
     19   PUBLIC solve_at_general_sed ! 
     20 
     21   ! Maximum number of iterations for each method 
     22   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20 
     23   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 
    1424 
    1525   !! * Module variables 
    1626   REAL(wp) :: & 
    17       salchl = 1./1.80655 , & ! conversion factor for salinity --> chlorinity (Wooster et al. 1969) 
    1827      calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]  
    1928 
    20    REAL(wp) :: &              ! coeff. for 1. dissoc. of carbonic acid (Millero, 1995)    
    21       c10 = 3670.7   , & 
    22       c11 = 62.008   , & 
    23       c12 = 9.7944   , & 
    24       c13 = 0.0118   , & 
    25       c14 = 0.000116 
    26  
    27    REAL(wp) :: &              ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)    
    28       c20 = 1394.7   , & 
    29       c21 = 4.777    , & 
    30       c22 = 0.       , & 
    31       c23 = 0.0184   , & 
    32       c24 = 0.000118 
    33  
    34    REAL(wp) :: &              ! coeff. for 1. dissoc. of boric acid (Dickson and Goyet, 1994) 
    35       cb0  = -8966.90, & 
    36       cb1  = -2890.53, & 
    37       cb2  = -77.942 , & 
    38       cb3  = 1.728   , & 
    39       cb4  = -0.0996 , & 
    40       cb5  = 148.0248, & 
    41       cb6  = 137.1942, & 
    42       cb7  = 1.62142 , & 
    43       cb8  = -24.4344, & 
    44       cb9  = -25.085 , & 
    45       cb10 = -0.2474 , & 
    46       cb11 = 0.053105 
    47  
    48    REAL(wp) :: &             ! borat constants 
    49       bor1 = 0.000232, & 
    50       bor2 = 1./10.811 
    51  
    52    REAL(wp) :: &             ! constants for calculate concentrations  
    53       st1  = 0.14     , &    ! for sulfate (Morris & Riley 1966) 
    54       st2  = 1./96.062, & 
    55       ks0  = 141.328  , & 
    56       ks1  = -4276.1  , & 
    57       ks2  = -23.093  , & 
    58       ks3  = -13856.  , & 
    59       ks4  = 324.57   , & 
    60       ks5  = -47.986  , & 
    61       ks6  = 35474.   , & 
    62       ks7  = -771.54  , & 
    63       ks8  = 114.723  , & 
    64       ks9  = -2698.   , & 
    65       ks10 = 1776.    , & 
    66       ks11 = 1.       , & 
    67       ks12 = -0.001005  
    68  
    69    REAL(wp) :: &             ! constants for calculate concentrations  
    70       ft1  = 0.000067   , &  ! fluorides (Dickson & Riley 1979 ) 
    71       ft2  = 1./18.9984 , & 
    72       kf0  = -12.641    , & 
    73       kf1  = 1590.2     , & 
    74       kf2  = 1.525      , & 
    75       kf3  = 1.0        , & 
    76       kf4  =-0.001005 
    77  
    78    REAL(wp) :: &             ! coeff. for dissoc. of water (Dickson and Riley, 1979 ) 
    79       cw0 = -13847.26  , & 
    80       cw1 = 148.9802   , & 
    81       cw2 = -23.6521   , & 
    82       cw3 = 118.67     , & 
    83       cw4 = -5.977     , & 
    84       cw5 = 1.0495     , & 
    85       cw6 = -0.01615 
    86   
    87    REAL(wp) :: &             ! coeff. for dissoc. of phosphate (Millero (1974) 
    88       cp10 = 115.54    , & 
    89       cp11 = -4576.752 , & 
    90       cp12 = -18.453   , & 
    91       cp13 = -106.736  , & 
    92       cp14 = 0.69171   , & 
    93       cp15 = -0.65643  , & 
    94       cp16 = -0.01844  , & 
    95       cp20 = 172.1033  , & 
    96       cp21 = -8814.715 , & 
    97       cp22 = -27.927   , & 
    98       cp23 = -160.340  , & 
    99       cp24 = 1.3566    , & 
    100       cp25 = 0.37335   , & 
    101       cp26 = -0.05778  , & 
    102       cp30 = -18.126   , & 
    103       cp31 = -3070.75  , & 
    104       cp32 = 17.27039  , & 
    105       cp33 = 2.81197   , & 
    106       cp34 = -44.99486 , & 
    107       cp35 = -0.09984 
    108  
    109    REAL(wp) :: &             ! coeff. for dissoc. of silicates (Millero (1974)   
    110       cs10 = 117.40    , &   
    111       cs11 = -8904.2   , &  
    112       cs12 = -19.334   , &  
    113       cs13 = -458.79   , &  
    114       cs14 = 3.5913    , &  
    115       cs15 = 188.74    , &  
    116       cs16 = -1.5998   , &  
    117       cs17 = -12.1652  , &  
    118       cs18 = 0.07871   , &  
    119       cs19 = 0.        , &  
    120       cs20 = 1.        , &  
    121       cs21 = -0.001005 
    122  
    123    REAL(wp) :: &            ! coeff. for apparent solubility equilibrium  
    124       ! akcc1 = -34.452 , & ! of calcite (Ingle,1975,  1800, eq. 6) 
    125       ! akcc2 = -39.866 , &   
    126       ! akcc3 = 110.21  , &   
    127       ! akcc4 = -7.5752E-6   
    128       akcc1 = -171.9065 , &    ! Millero et al. 1995 from Mucci 1983 
    129       akcc2 = -0.077993 , &   
    130       akcc3 = 2839.319  , &   
    131       akcc4 = 71.595    , &   
    132       akcc5 = -0.77712  , &   
    133       akcc6 = 0.0028426 , &   
    134       akcc7 = 178.34    , &   
    135       akcc8 = -0.07711  , &   
    136       akcc9 = 0.0041249 
    137  
    138    REAL(wp) :: &                 ! universal gas constant 
    139       rgas = 83.145              ! bar.cm3/(mol.K) or R=8.31451 J/(g.mol.K) 
    140  
    141  
    142    ! coeff. for seawater pressure correction (millero 95)       
    143    REAL(wp), DIMENSION(8)  :: &  
    144       devk1, devk2, devk3, devk4, devk5 
    145  
    146    DATA devk1/ -25.5   , -15.82    , -29.48  , -25.60    , -48.76    , -14.51   , -23.12   , -26.57   /    
    147    DATA devk2/ 0.1271  , -0.0219   , 0.1622  , 0.2324    , -0.5304   , 0.1211   ,  0.1758  , 0.2020   /    
    148    DATA devk3/ 0.      , 0.        , 2.608E-3, -3.6246E-3, 0.        , -0.321E-3, -2.647E-3, -3.042E-3/    
    149    DATA devk4/-3.08E-3 , 1.13E-3   , -2.84E-3, -5.13E-3  , -11.76E-3 , -2.67E-3 , -5.15E-3 , -4.08E-3 /    
    150    DATA devk5/0.0877E-3, -0.1475E-3, 0.      , 0.0794E-3 , -0.3692E-3, 0.0427E-3,  0.09E-3 , 0.0714E-3/ 
    151  
     29   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants 
    15230 
    15331   ! coeff. for density of sea water (Millero & Poisson 1981)  
     
    16240   REAL(wp), DIMENSION(6)  :: Ddsw                     
    16341   DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/ 
     42 
     43  REAL(wp) :: devk10  = -25.5 
     44   REAL(wp) :: devk11  = -15.82 
     45   REAL(wp) :: devk12  = -29.48 
     46   REAL(wp) :: devk13  = -20.02 
     47   REAL(wp) :: devk14  = -18.03 
     48   REAL(wp) :: devk15  = -9.78 
     49   REAL(wp) :: devk16  = -48.76 
     50   REAL(wp) :: devk17  = -14.51 
     51   REAL(wp) :: devk18  = -23.12 
     52   REAL(wp) :: devk19  = -26.57 
     53   REAL(wp) :: devk110  = -29.48 
     54   ! 
     55   REAL(wp) :: devk20  = 0.1271 
     56   REAL(wp) :: devk21  = -0.0219 
     57   REAL(wp) :: devk22  = 0.1622 
     58   REAL(wp) :: devk23  = 0.1119 
     59   REAL(wp) :: devk24  = 0.0466 
     60   REAL(wp) :: devk25  = -0.0090 
     61   REAL(wp) :: devk26  = 0.5304 
     62   REAL(wp) :: devk27  = 0.1211 
     63   REAL(wp) :: devk28  = 0.1758 
     64   REAL(wp) :: devk29  = 0.2020 
     65   REAL(wp) :: devk210  = 0.1622 
     66   ! 
     67   REAL(wp) :: devk30  = 0. 
     68   REAL(wp) :: devk31  = 0. 
     69   REAL(wp) :: devk32  = 2.608E-3 
     70   REAL(wp) :: devk33  = -1.409e-3 
     71   REAL(wp) :: devk34  = 0.316e-3 
     72   REAL(wp) :: devk35  = -0.942e-3 
     73   REAL(wp) :: devk36  = 0. 
     74   REAL(wp) :: devk37  = -0.321e-3 
     75   REAL(wp) :: devk38  = -2.647e-3 
     76   REAL(wp) :: devk39  = -3.042e-3 
     77   REAL(wp) :: devk310  = -2.6080e-3 
     78   ! 
     79   REAL(wp) :: devk40  = -3.08E-3 
     80   REAL(wp) :: devk41  = 1.13E-3 
     81   REAL(wp) :: devk42  = -2.84E-3 
     82   REAL(wp) :: devk43  = -5.13E-3 
     83   REAL(wp) :: devk44  = -4.53e-3 
     84   REAL(wp) :: devk45  = -3.91e-3 
     85   REAL(wp) :: devk46  = -11.76e-3 
     86   REAL(wp) :: devk47  = -2.67e-3 
     87   REAL(wp) :: devk48  = -5.15e-3 
     88   REAL(wp) :: devk49  = -4.08e-3 
     89   REAL(wp) :: devk410  = -2.84e-3 
     90   ! 
     91   REAL(wp) :: devk50  = 0.0877E-3 
     92   REAL(wp) :: devk51  = -0.1475E-3 
     93   REAL(wp) :: devk52  = 0. 
     94   REAL(wp) :: devk53  = 0.0794E-3 
     95   REAL(wp) :: devk54  = 0.09e-3 
     96   REAL(wp) :: devk55  = 0.054e-3 
     97   REAL(wp) :: devk56  = 0.3692E-3 
     98   REAL(wp) :: devk57  = 0.0427e-3 
     99   REAL(wp) :: devk58  = 0.09e-3 
     100   REAL(wp) :: devk59  = 0.0714e-3 
     101   REAL(wp) :: devk510  = 0.0 
    164102 
    165103   !! $Id$ 
     
    179117      INTEGER, INTENT(in) :: kt                     ! time step 
    180118 
    181 #if ! defined key_sed_off 
    182119      INTEGER  :: ji, jj, ikt 
    183       REAL(wp) :: ztkel, ztc, ztc2, zpres, ztr  
    184       REAL(wp) :: zsal, zsal2, zsqrt, zsal15   
    185       REAL(wp) :: zis, zis2, zisqrt           
     120      REAL(wp) :: ztc, ztc2 
     121      REAL(wp) :: zsal, zsal15   
    186122      REAL(wp) :: zdens0, zaw, zbw, zcw     
    187       REAL(wp) :: zbuf1, zbuf2  
    188       REAL(wp) :: zcpexp, zcpexp2 
    189       REAL(wp) :: zck1p, zck2p, zck3p, zcksi    
    190       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zchem_data 
    191  
    192 #endif 
     123      REAL(wp), DIMENSION(jpi,jpj,15) ::   zchem_data 
    193124      !!---------------------------------------------------------------------- 
    194125 
    195       IF( MOD( kt - 1, nfreq ) /= 0 ) RETURN 
    196  
    197       WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt 
    198       WRITE(numsed,*) ' ' 
    199  
    200  
    201 #if defined key_sed_off 
    202       CALL sed_chem_off 
    203 #else 
     126 
     127      IF( ln_timing )  CALL timing_start('sed_chem') 
     128 
     129      IF (lwp) WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt 
     130      IF (lwp) WRITE(numsed,*) ' ' 
     131 
    204132      ! reading variables 
    205       ALLOCATE( zchem_data(jpi,jpj,6) )   ;   zchem_data(:,:,:) = 0. 
    206  
    207       DO jj = 1,jpj 
    208          DO ji = 1, jpi 
    209             ikt = mbkt(ji,jj)  
    210             IF ( tmask(ji,jj,ikt) == 1 ) THEN 
    211                zchem_data(ji,jj,1) = ak13 (ji,jj,ikt) 
    212                zchem_data(ji,jj,2) = ak23 (ji,jj,ikt) 
    213                zchem_data(ji,jj,3) = akb3 (ji,jj,ikt) 
    214                zchem_data(ji,jj,4) = akw3 (ji,jj,ikt) 
    215                zchem_data(ji,jj,5) = aksp (ji,jj,ikt) 
    216                zchem_data(ji,jj,6) = borat(ji,jj,ikt) 
    217             ENDIF 
     133      zchem_data(:,:,:) = rtrn 
     134 
     135      IF (ln_sediment_offline) THEN 
     136         CALL sed_chem_cst 
     137      ELSE 
     138         DO jj = 1,jpj 
     139            DO ji = 1, jpi 
     140               ikt = mbkt(ji,jj)  
     141               IF ( tmask(ji,jj,ikt) == 1 ) THEN 
     142                  zchem_data(ji,jj,1) = ak13  (ji,jj,ikt) 
     143                  zchem_data(ji,jj,2) = ak23  (ji,jj,ikt) 
     144                  zchem_data(ji,jj,3) = akb3  (ji,jj,ikt) 
     145                  zchem_data(ji,jj,4) = akw3  (ji,jj,ikt) 
     146                  zchem_data(ji,jj,5) = aksp  (ji,jj,ikt) 
     147                  zchem_data(ji,jj,6) = borat (ji,jj,ikt) 
     148                  zchem_data(ji,jj,7) = ak1p3 (ji,jj,ikt) 
     149                  zchem_data(ji,jj,8) = ak2p3 (ji,jj,ikt) 
     150                  zchem_data(ji,jj,9) = ak3p3 (ji,jj,ikt) 
     151                  zchem_data(ji,jj,10)= aksi3 (ji,jj,ikt) 
     152                  zchem_data(ji,jj,11)= sio3eq(ji,jj,ikt) 
     153                  zchem_data(ji,jj,12)= aks3  (ji,jj,ikt) 
     154                  zchem_data(ji,jj,13)= akf3  (ji,jj,ikt) 
     155                  zchem_data(ji,jj,14)= sulfat(ji,jj,ikt) 
     156                  zchem_data(ji,jj,15)= fluorid(ji,jj,ikt) 
     157               ENDIF 
     158            ENDDO 
    218159         ENDDO 
    219       ENDDO 
    220  
    221       CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1), iarroce(1:jpoce) ) 
    222       CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2), iarroce(1:jpoce) ) 
    223       CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3), iarroce(1:jpoce) ) 
    224       CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4), iarroce(1:jpoce) ) 
    225       CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5), iarroce(1:jpoce) ) 
    226       CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6), iarroce(1:jpoce) ) 
    227  
    228       DEALLOCATE( zchem_data ) 
     160 
     161         CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1) , iarroce(1:jpoce) ) 
     162         CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2) , iarroce(1:jpoce) ) 
     163         CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3) , iarroce(1:jpoce) ) 
     164         CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4) , iarroce(1:jpoce) ) 
     165         CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5) , iarroce(1:jpoce) ) 
     166         CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6) , iarroce(1:jpoce) ) 
     167         CALL pack_arr ( jpoce, ak1ps (1:jpoce), zchem_data(1:jpi,1:jpj,7) , iarroce(1:jpoce) ) 
     168         CALL pack_arr ( jpoce, ak2ps (1:jpoce), zchem_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) ) 
     169         CALL pack_arr ( jpoce, ak3ps (1:jpoce), zchem_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) ) 
     170         CALL pack_arr ( jpoce, aksis (1:jpoce), zchem_data(1:jpi,1:jpj,10), iarroce(1:jpoce) ) 
     171         CALL pack_arr ( jpoce, sieqs (1:jpoce), zchem_data(1:jpi,1:jpj,11), iarroce(1:jpoce) ) 
     172         CALL pack_arr ( jpoce, aks3s (1:jpoce), zchem_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 
     173         CALL pack_arr ( jpoce, akf3s (1:jpoce), zchem_data(1:jpi,1:jpj,13), iarroce(1:jpoce) ) 
     174         CALL pack_arr ( jpoce, sulfats(1:jpoce), zchem_data(1:jpi,1:jpj,14), iarroce(1:jpoce) ) 
     175         CALL pack_arr ( jpoce, fluorids(1:jpoce), zchem_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 
     176      ENDIF 
    229177 
    230178      DO ji = 1, jpoce 
    231          ztkel   = temp(ji) + 273.16 
    232179         ztc     = temp(ji) 
    233180         ztc2    = ztc * ztc 
    234          zpres   = press(ji) 
    235181         ! zqtt    = ztkel * 0.01 
    236182         zsal    = salt(ji) 
    237          zsal2   = zsal * zsal  
    238          zsqrt   = SQRT( zsal ) 
    239          zsal15  = zsqrt * zsal 
    240          zlogt   = LOG( ztkel ) 
    241          ztr     = 1./ ztkel 
    242          ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet) 
    243          zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal ) 
    244          zis2    = zis * zis 
    245          zisqrt  = SQRT( zis ) 
     183         zsal15  = SQRT( zsal ) * zsal 
    246184 
    247185         ! Density of Sea Water - F(temp,sal) [kg/m3] 
     
    256194         densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l] 
    257195 
    258          ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)  
    259          ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982)  
    260          ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN 
    261          ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.)  
    262          ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP  
    263          ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286  
    264          ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285) 
    265          !----------------------------------------------------------------------------------------- 
    266          zcpexp  = zpres / ( rgas*ztkel ) 
    267          zcpexp2 = zpres * zcpexp 
    268  
    269          ! For Phodphates (phosphoric acid) (DOE 1994) 
    270          !---------------------------------------------- 
    271          zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt & 
    272             &      + ( cp15*ztr + cp16 ) * zsal 
    273          zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt & 
    274             &      + ( cp25*ztr + cp26 ) * zsal 
    275          zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt & 
    276             &      + ( cp34*ztr + cp35 ) * zsal 
    277  
    278          ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP) 
    279          !-------------------------------------------------------- 
    280          zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt & 
    281             &      + ( cs15*ztr + cs16 ) * zis & 
    282             &      + ( cs17*ztr + cs18 ) * zis2 & 
    283             &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal ) 
    284  
    285  
    286          !K1, K2 of carbonic acid, KB of boric acid, KW (H2O) 
    287          !--------------------------------------------------- 
    288          zak1p  = EXP ( zck1p  ) 
    289          zak2p  = EXP ( zck2p  ) 
    290          zak3p  = EXP ( zck3p  ) 
    291          zaksil = EXP ( zcksi  ) 
    292  
    293          zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 ) 
    294          zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc ) 
    295          aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    296  
    297          zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 ) 
    298          zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc ) 
    299          ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    300   
    301          zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 ) 
    302          zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc ) 
    303          ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    304  
    305          zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 ) 
    306          zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc ) 
    307          ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    308  
    309196         ak12s  (ji) = ak1s (ji) * ak2s (ji) 
    310197         ak12ps (ji) = ak1ps(ji) * ak2ps(ji) 
     
    313200         calcon2(ji) = 0.01028 * ( salt(ji) / 35. ) * densSW(ji) 
    314201      ENDDO 
    315  
    316202        
    317 #endif 
     203      IF( ln_timing )  CALL timing_stop('sed_chem') 
    318204 
    319205   END SUBROUTINE sed_chem 
    320206 
    321 #if defined key_sed_off 
    322  
    323    SUBROUTINE sed_chem_off 
    324       !!---------------------------------------------------------------------- 
    325       !!                   ***  ROUTINE sed_chem_off  *** 
    326       !!                     
    327       !! ** Purpose :   compute chemical constants 
     207   SUBROUTINE ahini_for_at_sed(p_hini) 
     208      !!--------------------------------------------------------------------- 
     209      !!                     ***  ROUTINE ahini_for_at  *** 
    328210      !! 
    329       !!   History : 
    330       !!        !  04-10  (N. Emprin, M. Gehlen )  Original code 
    331       !!        !  06-04  (C. Ethe)  Re-organization  
    332       !!----------------------------------------------------------------------   
    333       !! * Local declarations 
    334       INTEGER :: ji 
    335  
    336       REAL(wp) ::  ztkel, ztc, ztc2, zpres, ztr  
    337       REAL(wp) ::  zsal, zsal2, zsqrt, zsal15    
    338       REAL(wp) ::  zis, zis2, zisqrt             
    339       REAL(wp) ::  zdens0, zaw, zbw, zcw         
    340       REAL(wp) ::  zchl, zst, zft, zbuf1, zbuf2  
    341       REAL(wp) ::  zcpexp, zcpexp2               
    342       REAL(wp) ::  zckb, zck1, zck2, zckw        
    343       REAL(wp) ::  zck1p, zck2p, zck3p, zcksi    
    344       REAL(wp) ::  zak1, zak2, zakb, zakw               
    345       REAL(wp) ::  zaksp0, zksp, zks, zkf  
    346  
     211      !! Subroutine returns the root for the 2nd order approximation of the 
     212      !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic  
     213      !! polynomial) around the local minimum, if it exists. 
     214      !! Returns * 1E-03_wp if p_alkcb <= 0 
     215      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 
     216      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 
     217      !!                    and the 2nd order approximation does not have  
     218      !!                    a solution 
     219      !!--------------------------------------------------------------------- 
     220      REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  ::  p_hini 
     221      INTEGER  ::   ji, jk 
     222      REAL(wp)  ::  zca1, zba1 
     223      REAL(wp)  ::  zd, zsqrtd, zhmin 
     224      REAL(wp)  ::  za2, za1, za0 
     225      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb 
     226 
     227      IF( ln_timing )  CALL timing_start('ahini_for_at_sed') 
     228      ! 
     229      DO jk = 1, jpksed 
     230         DO ji = 1, jpoce 
     231            p_alkcb  = pwcp(ji,jk,jwalk) / densSW(ji) 
     232            p_dictot = pwcp(ji,jk,jwdic) / densSW(ji)  
     233            p_bortot = borats(ji) / densSW(ji) 
     234            IF (p_alkcb <= 0.) THEN 
     235                p_hini(ji,jk) = 1.e-3 
     236            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 
     237                p_hini(ji,jk) = 1.e-10_wp 
     238            ELSE 
     239                zca1 = p_dictot/( p_alkcb + rtrn ) 
     240                zba1 = p_bortot/ (p_alkcb + rtrn ) 
     241           ! Coefficients of the cubic polynomial 
     242                za2 = aKbs(ji)*(1. - zba1) + ak1s(ji)*(1.-zca1) 
     243                za1 = ak1s(ji)*akbs(ji)*(1. - zba1 - zca1)    & 
     244                &     + ak1s(ji)*ak2s(ji)*(1. - (zca1+zca1)) 
     245                za0 = ak1s(ji)*ak2s(ji)*akbs(ji)*(1. - zba1 - (zca1+zca1)) 
     246                                        ! Taylor expansion around the minimum 
     247                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation 
     248                                        ! for the minimum close to the root 
     249 
     250                IF(zd > 0.) THEN        ! If the discriminant is positive 
     251                  zsqrtd = SQRT(zd) 
     252                  IF(za2 < 0) THEN 
     253                    zhmin = (-za2 + zsqrtd)/3. 
     254                  ELSE 
     255                    zhmin = -za1/(za2 + zsqrtd) 
     256                  ENDIF 
     257                  p_hini(ji,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 
     258                ELSE 
     259                  p_hini(ji,jk) = 1.e-7 
     260                ENDIF 
     261             ! 
     262            ENDIF 
     263         END DO 
     264      END DO 
     265      ! 
     266      IF( ln_timing )  CALL timing_stop('ahini_for_at_sed') 
     267      ! 
     268   END SUBROUTINE ahini_for_at_sed 
     269 
     270   !=============================================================================== 
     271   SUBROUTINE anw_infsup_sed( p_alknw_inf, p_alknw_sup ) 
     272 
     273   ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 
     274   ! contributions to total alkalinity (the infimum and the supremum), i.e 
     275   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 
     276 
     277   ! Argument variables 
     278   INTEGER :: jk 
     279   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_inf 
     280   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_sup 
     281 
     282   DO jk = 1, jpksed 
     283      p_alknw_inf(:,jk) =  -pwcp(:,jk,jwpo4) / densSW(:) 
     284      p_alknw_sup(:,jk) =   (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil)     & 
     285   &                        + borats(:) ) / densSW(:) 
     286   END DO 
     287 
     288   END SUBROUTINE anw_infsup_sed 
     289 
     290 
     291   SUBROUTINE solve_at_general_sed( p_hini, zhi ) 
     292 
     293   ! Universal pH solver that converges from any given initial value, 
     294   ! determines upper an lower bounds for the solution if required 
     295 
     296   ! Argument variables 
     297   !-------------------- 
     298   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(IN)   :: p_hini 
     299   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  :: zhi 
     300 
     301   ! Local variables 
     302   !----------------- 
     303   INTEGER   ::  ji, jk, jn 
     304   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor 
     305   REAL(wp)  ::  zdelta, zh_delta 
     306   REAL(wp)  ::  zeqn, zdeqndh, zalka 
     307   REAL(wp)  ::  aphscale 
     308   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 
     309   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 
     310   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 
     311   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 
     312   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 
     313   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 
     314   REAL(wp)  ::  zalk_wat, zdalk_wat 
     315   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 
     316   LOGICAL   ::  l_exitnow 
     317   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 
     318   REAL(wp), DIMENSION(jpoce,jpksed) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 
     319 
     320   IF( ln_timing )  CALL timing_start('solve_at_general_sed') 
     321      !  Allocate temporary workspace 
     322   CALL anw_infsup_sed( zalknw_inf, zalknw_sup ) 
     323 
     324   rmask(:,:) = 1.0 
     325   zhi(:,:)   = 0. 
     326 
     327   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 
     328   DO jk = 1, jpksed 
     329      DO ji = 1, jpoce 
     330         IF (rmask(ji,jk) == 1.) THEN 
     331            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji) 
     332            aphscale = 1. + sulfats(ji)/aks3s(ji) 
     333            zh_ini = p_hini(ji,jk) 
     334 
     335            zdelta = (p_alktot-zalknw_inf(ji,jk))**2 + 4.*akws(ji) / aphscale 
     336 
     337            IF(p_alktot >= zalknw_inf(ji,jk)) THEN 
     338               zh_min(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_inf(ji,jk) + SQRT(zdelta) ) 
     339            ELSE 
     340               zh_min(ji,jk) = aphscale * (-(p_alktot-zalknw_inf(ji,jk)) + SQRT(zdelta) ) / 2. 
     341            ENDIF 
     342 
     343            zdelta = (p_alktot-zalknw_sup(ji,jk))**2 + 4.*akws(ji) / aphscale 
     344 
     345            IF(p_alktot <= zalknw_sup(ji,jk)) THEN 
     346               zh_max(ji,jk) = aphscale * (-(p_alktot-zalknw_sup(ji,jk)) + SQRT(zdelta) ) / 2. 
     347            ELSE 
     348               zh_max(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_sup(ji,jk) + SQRT(zdelta) ) 
     349            ENDIF 
     350 
     351            zhi(ji,jk) = MAX(MIN(zh_max(ji,jk), zh_ini), zh_min(ji,jk)) 
     352         ENDIF 
     353      END DO 
     354   END DO 
     355 
     356   zeqn_absmin(:,:) = HUGE(1._wp) 
     357 
     358   DO jn = 1, jp_maxniter_atgen 
     359   DO jk = 1, jpksed 
     360      DO ji = 1, jpoce 
     361         IF (rmask(ji,jk) == 1.) THEN 
     362 
     363            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji) 
     364            zdic  = pwcp(ji,jk,jwdic) / densSW(ji) 
     365            zbot  = borats(ji) / densSW(ji) 
     366            zpt =  pwcp(ji,jk,jwpo4) / densSW(ji) 
     367            zsit = pwcp(ji,jk,jwsil) / densSW(ji) 
     368            zst = sulfats(ji) 
     369            zft = fluorids(ji) 
     370            aphscale = 1. + sulfats(ji)/aks3s(ji) 
     371            zh = zhi(ji,jk) 
     372            zh_prev = zh 
     373 
     374            ! H2CO3 - HCO3 - CO3 : n=2, m=0 
     375            znumer_dic = 2.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji) 
     376            zdenom_dic = ak1s(ji)*ak2s(ji) + zh*(ak1s(ji) + zh) 
     377            zalk_dic   = zdic * (znumer_dic/zdenom_dic) 
     378            zdnumer_dic = ak1s(ji)*ak1s(ji)*ak2s(ji) + zh     & 
     379                          *(4.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji)) 
     380            zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2) 
     381 
     382 
     383            ! B(OH)3 - B(OH)4 : n=1, m=0 
     384            znumer_bor = akbs(ji) 
     385            zdenom_bor = akbs(ji) + zh 
     386            zalk_bor   = zbot * (znumer_bor/zdenom_bor) 
     387            zdnumer_bor = akbs(ji) 
     388            zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2) 
     389 
     390 
     391            ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 
     392            znumer_po4 = 3.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  & 
     393            &            + zh*(2.*ak1ps(ji)*ak2ps(ji) + zh* ak1ps(ji)) 
     394            zdenom_po4 = ak1ps(ji)*ak2ps(ji)*ak3ps(ji)     & 
     395            &            + zh*( ak1ps(ji)*ak2ps(ji) + zh*(ak1ps(ji) + zh)) 
     396            zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 
     397            zdnumer_po4 = ak1ps(ji)*ak2ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  & 
     398            &             + zh*(4.*ak1ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)         & 
     399            &             + zh*(9.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)                         & 
     400            &             + ak1ps(ji)*ak1ps(ji)*ak2ps(ji)                                & 
     401            &             + zh*(4.*ak1ps(ji)*ak2ps(ji) + zh * ak1ps(ji) ) ) ) 
     402            zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2) 
     403 
     404            ! H4SiO4 - H3SiO4 : n=1, m=0 
     405            znumer_sil = aksis(ji) 
     406            zdenom_sil = aksis(ji) + zh 
     407            zalk_sil   = zsit * (znumer_sil/zdenom_sil) 
     408            zdnumer_sil = aksis(ji) 
     409            zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2) 
     410 
     411            ! HSO4 - SO4 : n=1, m=1 
     412            aphscale = 1.0 + zst/aks3s(ji) 
     413            znumer_so4 = aks3s(ji) * aphscale 
     414            zdenom_so4 = aks3s(ji) * aphscale + zh 
     415            zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.) 
     416            zdnumer_so4 = aks3s(ji) * aphscale 
     417            zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2) 
     418 
     419            ! HF - F : n=1, m=1 
     420            znumer_flu =  akf3s(ji) 
     421            zdenom_flu =  akf3s(ji) + zh 
     422            zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.) 
     423            zdnumer_flu = akf3s(ji) 
     424            zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2) 
     425 
     426            ! H2O - OH 
     427            zalk_wat   = akws(ji)/zh - zh/aphscale 
     428            zdalk_wat  = -akws(ji)/zh**2 - 1./aphscale 
     429 
     430            ! CALCULATE [ALK]([CO3--], [HCO3-]) 
     431            zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   & 
     432            &      + zalk_so4 + zalk_flu                       & 
     433            &      + zalk_wat - p_alktot 
     434 
     435            zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   & 
     436            &       + zalk_so4 + zalk_flu + zalk_wat) 
     437 
     438            zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 
     439            &         + zdalk_so4 + zdalk_flu + zdalk_wat 
     440 
     441            ! Adapt bracketing interval 
     442            IF(zeqn > 0._wp) THEN 
     443               zh_min(ji,jk) = zh_prev 
     444            ELSEIF(zeqn < 0._wp) THEN 
     445               zh_max(ji,jk) = zh_prev 
     446            ENDIF 
     447 
     448            IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jk)) THEN 
     449            ! if the function evaluation at the current point is 
     450            ! not decreasing faster than with a bisection step (at least linearly) 
     451            ! in absolute value take one bisection step on [ph_min, ph_max] 
     452            ! ph_new = (ph_min + ph_max)/2d0 
     453            ! 
     454            ! In terms of [H]_new: 
     455            ! [H]_new = 10**(-ph_new) 
     456            !         = 10**(-(ph_min + ph_max)/2d0) 
     457            !         = SQRT(10**(-(ph_min + phmax))) 
     458            !         = SQRT(zh_max * zh_min) 
     459               zh = SQRT(zh_max(ji,jk) * zh_min(ji,jk)) 
     460               zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     461            ELSE 
     462            ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 
     463            !           = -zdeqndh * LOG(10) * [H] 
     464            ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 
     465            ! 
     466            ! pH_new = pH_old + \deltapH 
     467            ! 
     468            ! [H]_new = 10**(-pH_new) 
     469            !         = 10**(-pH_old - \Delta pH) 
     470            !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 
     471            !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 
     472            !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 
     473 
     474               zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 
     475 
     476               IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 
     477                  zh          = zh_prev*EXP(zh_lnfactor) 
     478               ELSE 
     479                  zh_delta    = zh_lnfactor*zh_prev 
     480                  zh          = zh_prev + zh_delta 
     481               ENDIF 
     482 
     483               IF( zh < zh_min(ji,jk) ) THEN 
     484               ! if [H]_new < [H]_min 
     485               ! i.e., if ph_new > ph_max then 
     486               ! take one bisection step on [ph_prev, ph_max] 
     487               ! ph_new = (ph_prev + ph_max)/2d0 
     488               ! In terms of [H]_new: 
     489               ! [H]_new = 10**(-ph_new) 
     490               !         = 10**(-(ph_prev + ph_max)/2d0) 
     491               !         = SQRT(10**(-(ph_prev + phmax))) 
     492               !         = SQRT([H]_old*10**(-ph_max)) 
     493               !         = SQRT([H]_old * zh_min) 
     494                  zh                = SQRT(zh_prev * zh_min(ji,jk)) 
     495                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     496               ENDIF 
     497 
     498               IF( zh > zh_max(ji,jk) ) THEN 
     499               ! if [H]_new > [H]_max 
     500               ! i.e., if ph_new < ph_min, then 
     501               ! take one bisection step on [ph_min, ph_prev] 
     502               ! ph_new = (ph_prev + ph_min)/2d0 
     503               ! In terms of [H]_new: 
     504               ! [H]_new = 10**(-ph_new) 
     505               !         = 10**(-(ph_prev + ph_min)/2d0) 
     506               !         = SQRT(10**(-(ph_prev + ph_min))) 
     507               !         = SQRT([H]_old*10**(-ph_min)) 
     508               !         = SQRT([H]_old * zhmax) 
     509                  zh                = SQRT(zh_prev * zh_max(ji,jk)) 
     510                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     511               ENDIF 
     512            ENDIF 
     513 
     514            zeqn_absmin(ji,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jk)) 
     515 
     516            ! Stop iterations once |\delta{[H]}/[H]| < rdel 
     517            ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 
     518            ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 
     519            ! Alternatively: 
     520            ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 
     521            !             ~ 1/LOG(10) * |\Delta [H]|/[H] 
     522            !             < 1/LOG(10) * rdel 
     523 
     524            ! Hence |zeqn/(zdeqndh*zh)| < rdel 
     525 
     526            ! rdel <-- pp_rdel_ah_target 
     527            l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 
     528 
     529            IF(l_exitnow) THEN 
     530               rmask(ji,jk) = 0. 
     531            ENDIF 
     532 
     533            zhi(ji,jk) =  zh 
     534 
     535            IF(jn >= jp_maxniter_atgen) THEN 
     536               zhi(ji,jk) = -1._wp 
     537            ENDIF 
     538 
     539         ENDIF 
     540      END DO 
     541   END DO 
     542   END DO 
     543   ! 
     544   IF( ln_timing )  CALL timing_stop('solve_at_general_sed') 
     545 
     546   END SUBROUTINE solve_at_general_sed 
     547 
     548   SUBROUTINE sed_chem_cst 
     549      !!--------------------------------------------------------------------- 
     550      !!                     ***  ROUTINE sed_chem_cst  *** 
     551      !! 
     552      !! ** Purpose :   Sea water chemistry computed following MOCSY protocol 
     553      !!                Computation is done at the bottom of the ocean only 
     554      !! 
     555      !! ** Method  : - ... 
     556      !!--------------------------------------------------------------------- 
     557      INTEGER  ::   ji 
     558      REAL(wp), DIMENSION(jpoce) :: saltprac, temps 
     559      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2 
     560      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 
     561      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2 
     562      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat 
     563      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2 
     564      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
     565      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 
     566      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1 
     567      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total 
     568      !!--------------------------------------------------------------------- 
     569      ! 
     570      IF( ln_timing )   CALL timing_start('sed_chem_cst') 
     571      ! 
     572      ! Computation of chemical constants require practical salinity 
     573      ! Thus, when TEOS08 is used, absolute salinity is converted to 
     574      ! practical salinity 
     575      ! ------------------------------------------------------------- 
     576      IF (neos == -1) THEN 
     577         saltprac(:) = salt(:) * 35.0 / 35.16504 
     578      ELSE 
     579         saltprac(:) = temp(:) 
     580      ENDIF 
     581 
     582      ! 
     583      ! Computations of chemical constants require in situ temperature 
     584      ! Here a quite simple formulation is used to convert 
     585      ! potential temperature to in situ temperature. The errors is less than 
     586      ! 0.04°C relative to an exact computation 
     587      ! --------------------------------------------------------------------- 
     588         DO ji = 1, jpoce 
     589            zpres = zkbot(ji) / 1000. 
     590            za1 = 0.04 * ( 1.0 + 0.185 * temp(ji) + 0.035 * (saltprac(ji) - 35.0) ) 
     591            za2 = 0.0075 * ( 1.0 - temp(ji) / 30.0 ) 
     592            temps(ji) = temp(ji) - za1 * zpres + za2 * zpres**2 
     593         END DO 
    347594 
    348595      ! CHEMICAL CONSTANTS - DEEP OCEAN 
    349       !------------------------------------- 
    350       ! [chem constants]=mol/kg solution (or (mol/kg sol)2 for akws and aksp) 
    351  
     596      ! ------------------------------- 
    352597      DO ji = 1, jpoce 
    353          ztkel   = temp(ji) + 273.16 
    354          ztc     = temp(ji) 
    355          ztc2    = ztc * ztc 
    356          zpres   = press(ji) 
    357          ! zqtt    = ztkel * 0.01 
    358          zsal    = salt(ji) 
    359          zsal2   = zsal * zsal  
    360          zsqrt   = SQRT( zsal ) 
     598         ! SET PRESSION ACCORDING TO SAUNDER (1980) 
     599         zc1 = 5.92E-3  
     600         zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*zkbot(ji)))) / 4.42E-6 
     601         zpres = zpres / 10.0 
     602 
     603         ! SET ABSOLUTE TEMPERATURE 
     604         ztkel   = temps(ji) + 273.15 
     605         zsal    = saltprac(ji) 
     606         zsqrt  = SQRT( zsal ) 
    361607         zsal15  = zsqrt * zsal 
    362          zlogt   = LOG( ztkel ) 
    363          ztr     = 1./ ztkel 
    364          ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet) 
    365          zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal ) 
    366          zis2    = zis * zis 
    367          zisqrt  = SQRT( zis ) 
    368  
    369  
    370          ! Density of Sea Water - F(temp,sal) [kg/m3] 
    371          zdens0 =  Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 & 
    372                   + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 & 
    373                   + Ddsw(6) * ztc * ztc2 * ztc2 
    374          zaw =  Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 & 
    375               + Adsw(5) * ztc2 * ztc2 
    376          zbw =  Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2 
    377          zcw =  Cdsw 
    378          densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal 
    379          densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l] 
    380  
    381  
    382          ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)  
    383          ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982)  
    384          ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN 
    385          ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.)  
    386          ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP  
    387          ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286  
    388          ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285) 
    389          !----------------------------------------------------------------------------------------- 
    390          zcpexp  = zpres / ( rgas*ztkel ) 
     608         zlogt  = LOG( ztkel ) 
     609         ztr    = 1. / ztkel 
     610         zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal ) 
     611         zis2   = zis * zis 
     612         zisqrt = SQRT( zis ) 
     613         ztc    = temps(ji) 
     614 
     615         ! CHLORINITY (WOOSTER ET AL., 1969) 
     616         zcl     = zsal / 1.80655 
     617 
     618         ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 
     619         zst     = 0.14 * zcl /96.062 
     620 
     621         ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 
     622         zft     = 0.000067 * zcl /18.9984 
     623 
     624         ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
     625         zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         & 
     626         &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 
     627         &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    & 
     628         &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         & 
     629         &         + LOG(1.0 - 0.001005 * zsal)) 
     630 
     631         ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
     632         zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   & 
     633         &         + LOG(1.0d0 - 0.001005d0*zsal)            & 
     634         &         + LOG(1.0d0 + zst/zcks)) 
     635 
     636         ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
     637         zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        & 
     638         &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         & 
     639         &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   & 
     640         &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
     641         &      * zlogt + 0.053105*zsqrt*ztkel 
     642 
     643         ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
     644         ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 
     645         zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  & 
     646                   - 0.011555*zsal + 0.0001152*zsal*zsal) 
     647         zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      & 
     648                   - 0.01781*zsal + 0.0001122*zsal*zsal) 
     649 
     650         ! PKW (H2O) (MILLERO, 1995) from composite data 
     651         zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    & 
     652                   - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 
     653 
     654         ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 
     655         zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   & 
     656         &          + (-106.736*ztr + 0.69171) * zsqrt       & 
     657         &          + (-0.65643*ztr - 0.01844) * zsal 
     658 
     659         zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  & 
     660         &          + (-160.340*ztr + 1.3566)*zsqrt          & 
     661         &          + (0.37335*ztr - 0.05778)*zsal 
     662 
     663         zck3p    = -3070.75*ztr - 18.126                    & 
     664         &          + (17.27039*ztr + 2.81197) * zsqrt       & 
     665         &          + (-44.99486*ztr - 0.09984) * zsal 
     666 
     667         ! CONSTANT FOR SILICATE, MILLERO (1995) 
     668         zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   & 
     669         &          + (-458.79*ztr + 3.5913) * zisqrt       & 
     670         &          + (188.74*ztr - 1.5998) * zis           & 
     671         &          + (-12.1652*ztr + 0.07871) * zis2       & 
     672         &          + LOG(1.0 - 0.001005*zsal) 
     673 
     674         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
     675         !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 
     676         zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   & 
     677         &         + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  & 
     678         &         - 0.07711*zsal + 0.0041249*zsal15 
     679 
     680         ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
     681         zak1    = 10**(zck1) * total2SWS 
     682         zak2    = 10**(zck2) * total2SWS 
     683         zakb    = EXP( zckb ) * total2SWS 
     684         zakw    = EXP( zckw ) 
     685         zaksp1  = 10**(zaksp0) 
     686         zak1p   = exp( zck1p ) 
     687         zak2p   = exp( zck2p ) 
     688         zak3p   = exp( zck3p ) 
     689         zaksi   = exp( zcksi ) 
     690         zckf    = zckf * total2SWS 
     691 
     692         ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 
     693         !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE 
     694         !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 
     695         !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN 
     696         !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 
     697         !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 
     698         !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 
     699         !        & GIESKES (1970), P. 1285-1286 (THE SMALL 
     700         !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 
     701         !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 
     702         zcpexp  = zpres / (rgas*ztkel) 
    391703         zcpexp2 = zpres * zcpexp 
    392704 
    393  
    394          ! chlorinity (WOOSTER ET AL., 1969) 
    395          !--------------------------------------- 
    396          zchl = zsal * salchl 
    397  
    398          ! total sulfate concentration [mol/kg soln] 
    399          ! -------------------------------------- 
    400          zst = st1 * zchl * st2 
    401  
    402          ! total fluoride concentration [mol/kg soln] 
    403          ! -------------------------------------- 
    404          zft = ft1 * zchl * ft2 
    405  
    406          ! dissociation constant for carbonate (Mehrback 74 - Dickson & Millero 87) 
    407          !--------------------------------------------------------------------------- 
    408          zck1 = c10*ztr - c11 + c12*zlogt - c13*zsal + c14*zsal2 
    409          zck2 = c20*ztr + c21 - c22*zlogt - c23*zsal + c24*zsal2 
    410  
    411          ! dissociation constant for sulfates (Dickson 1990) 
    412          !-------------------------------------------------- 
    413          zks = EXP(  ks0 + ks1*ztr + ks2*zlogt & 
    414             &    + ( ks3*ztr + ks4 + ks5*zlogt ) * zisqrt & 
    415             &    + ( ks6*ztr + ks7 + ks8*zlogt ) * zis    & 
    416             &    +   ks9*ztr*zis*zisqrt + ks10*ztr*zis2   & 
    417             &    +   LOG( ks11 + ks12*zsal ) ) 
    418  
    419          ! dissociation constant for fluorides (Dickson and Riley 79) 
    420          !-------------------------------------------------- 
    421          zkf = EXP( kf0 + kf1*ztr + kf2*zisqrt + LOG( kf3 + kf4*zsal ) ) 
    422  
    423          ! dissociation constant for borates (Doe 94) 
    424          !-------------------------------------------------- 
    425          zckb = (  cb0 + cb1*zsqrt + cb2*zsal + cb3*zsal15 + cb4*zsal2) * ztr & 
    426             &  + ( cb5 + cb6*zsqrt + cb7*zsal) & 
    427             &  + ( cb8 + cb9*zsqrt + cb10*zsal) * zlogt & 
    428             &  +   cb11*zsqrt*ztkel + LOG( ( 1. + zst/zks + zft/zkf ) / ( 1. + zst/zks ) )  
    429  
    430          ! PKW (H2O) (DICKSON AND RILEY, 1979) 
    431          !-------------------------------------- 
    432          zckw =   cw0*ztr + cw1 + cw2*zlogt & 
    433             & +( cw3*ztr + cw4 + cw5*zlogt )* zsqrt + cw6*zsal 
    434           
    435          ! For Phodphates (phosphoric acid) (DOE 1994) 
    436          !---------------------------------------------- 
    437          zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt & 
    438             &      + ( cp15*ztr + cp16 ) * zsal 
    439          zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt & 
    440             &      + ( cp25*ztr + cp26 ) * zsal 
    441          zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt & 
    442             &      + ( cp34*ztr + cp35 ) * zsal 
    443  
    444          ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP) 
    445          !-------------------------------------------------------- 
    446          zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt & 
    447             &      + ( cs15*ztr + cs16 ) * zis & 
    448             &      + ( cs17*ztr + cs18 ) * zis2 & 
    449             &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal ) 
    450  
    451          ! apparent solublity product K'SP of calcite in seawater 
    452          ! (S=27-43, T=2-25 DEG C) AT pres =0 (INGLE, 1975, EQ. 6) 
    453          ! prob: olivier a log = ln et C. Heize a LOG10(sal) 
    454          ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)+akcc4*tkel*tkel) 
    455          ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)+akcc4*tkel*tkel) 
    456          !-------------------------------------------------------------------- 
    457          zaksp0 = akcc1 + akcc2*ztkel + akcc3*ztr + akcc4 * LOG10(ztkel) & 
    458             &  + ( akcc5 + akcc6*ztkel+ akcc7*ztr ) * zsqrt & 
    459             &  +  akcc8*zsal + akcc9*zsal15 
    460  
    461          !K1, K2 of carbonic acid, KB of boric acid, KW (H2O) 
    462          !--------------------------------------------------- 
    463          zak1   = 10**( -zck1  ) 
    464          zak2   = 10**( -zck2  ) 
    465          zakb   = EXP ( zckb   )  
    466          zakw   = EXP ( zckw   ) 
    467          zksp   = 10**( zaksp0 ) 
    468  
    469  
    470  
    471          ! KB of boric acid, K1,K2 of carbonic acid pressure correction  
    472          ! after Culberson and  AND Pytkowicz (1968) (CF. BROECKER ET AL., 1982) Millero 95 
    473          !-------------------------------------------------------------------------------- 
    474          zbuf1       = - ( devk1(1) + devk2(1)*ztc + devk3(1)*ztc2 ) 
    475          zbuf2       = 0.5 * ( devk4(1) + devk5(1)*ztc ) 
    476          ak1s(ji)    = zak1 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    477  
    478          zbuf1       = -( devk1(2) + devk2(2)*ztc + devk3(2)*ztc2 ) 
    479          zbuf2       = 0.5 * ( devk4(2) + devk5(2)*ztc ) 
    480          ak2s(ji)    = zak2 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    481  
    482          zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 ) 
    483          zbuf2       = 0.5 * ( devk4(3) + devk5(3) * ztc ) 
    484          akbs(ji)    = zakb * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    485  
    486          zbuf1       = - ( devk1(4) + devk2(4)*ztc + devk3(4)*ztc2 ) 
    487          zbuf2       = 0.5 * ( devk4(4) + devk5(4)*ztc ) 
    488          akws(ji)    = zakw * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    489  
    490  
    491          ! APPARENT SOLUBILITY PRODUCT K''SP OF CALCITE (OR ARAGONITE) 
    492          ! AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970) 
    493          ! (P. 1285) AND BERNER (1976) 
    494          !----------------------------------------------------------------- 
    495          ! aksp(ji) = aksp0*exp(zcpexp*(devks-devkst*tc)) 
    496          ! or following Mucci 
    497          zbuf1      = - ( devk1(5) + devk2(5)*ztc + devk3(5)*ztc2 ) 
    498          zbuf2      = 0.5 *( devk4(5) + devk5(5)*ztc ) 
    499          aksps(ji)   = zksp * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    500  
    501          ! For Phodphates (phosphoric acid) (DOE 1994) 
    502          !---------------------------------------------- 
    503          zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt & 
    504             &      + ( cp15*ztr + cp16 ) * zsal 
    505          zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt & 
    506             &      + ( cp25*ztr + cp26 ) * zsal 
    507          zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt & 
    508             &      + ( cp34*ztr + cp35 ) * zsal 
    509  
    510          ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP) 
    511          !-------------------------------------------------------- 
    512          zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt & 
    513             &      + ( cs15*ztr + cs16 ) * zis & 
    514             &      + ( cs17*ztr + cs18 ) * zis2 & 
    515             &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal ) 
    516  
    517  
    518          !K1, K2 of carbonic acid, KB of boric acid, KW (H2O) 
    519          !--------------------------------------------------- 
    520          zak1p  = EXP ( zck1p  ) 
    521          zak2p  = EXP ( zck2p  ) 
    522          zak3p  = EXP ( zck3p  ) 
    523          zaksil = EXP ( zcksi  ) 
    524  
    525          zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 ) 
    526          zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc ) 
    527          aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    528  
    529          zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 ) 
    530          zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc ) 
    531          ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    532   
    533          zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 ) 
    534          zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc ) 
    535          ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    536  
    537          zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 ) 
    538          zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc ) 
    539          ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 
    540  
    541          ! total borat concentration. [mol/l] 
    542          ! or from Millero 1995 [mol/l] : borat(l) = 0.000416_8*(sal/35._8)*densSW(l) 
    543          ! -------------------------------------------------------------------------- 
    544          borats(ji) = bor1 * zchl * bor2 * densSW(ji) 
    545  
    546          ak12s  (ji) = ak1s (ji) * ak2s (ji) 
    547          ak12ps (ji) = ak1ps(ji) * ak2ps(ji) 
    548          ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji) 
    549  
    550          calcon2(ji) = 0.01028 * ( zsal / 35. ) * densSW(ji) 
    551  
    552       ENDDO 
    553  
    554    END SUBROUTINE sed_chem_off 
    555  
    556 #endif 
    557  
    558 #else 
    559    !!====================================================================== 
    560    !! MODULE sedchem  :   Dummy module  
    561    !!====================================================================== 
    562    !! $Id$ 
    563 CONTAINS 
    564    SUBROUTINE sed_chem( kt )         ! Empty routine 
    565       INTEGER, INTENT(in) :: kt 
    566       WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt 
    567    END SUBROUTINE sed_chem 
    568  
    569    !!====================================================================== 
    570  
    571 #endif 
     705         ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 
     706         !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 
     707         !        (CF. BROECKER ET AL., 1982) 
     708 
     709         zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 
     710         zbuf2  = 0.5 * ( devk40 + devk50 * ztc ) 
     711         ak1s(ji) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     712 
     713         zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 
     714         zbuf2  = 0.5 * ( devk41 + devk51 * ztc ) 
     715         ak2s(ji) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     716 
     717         zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 
     718         zbuf2  = 0.5 * ( devk42 + devk52 * ztc ) 
     719         akbs(ji) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     720 
     721         zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 
     722         zbuf2  = 0.5 * ( devk43 + devk53 * ztc ) 
     723         akws(ji) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     724 
     725         zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 
     726         zbuf2  = 0.5 * ( devk44 + devk54 * ztc ) 
     727         aks3s(ji) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     728 
     729         zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 
     730         zbuf2  = 0.5 * ( devk45 + devk55 * ztc ) 
     731         akf3s(ji) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     732 
     733         zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 
     734         zbuf2  = 0.5 * ( devk47 + devk57 * ztc ) 
     735         ak1ps(ji) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     736 
     737         zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 
     738         zbuf2  = 0.5 * ( devk48 + devk58 * ztc ) 
     739         ak2ps(ji) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     740 
     741         zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 
     742         zbuf2  = 0.5 * ( devk410 + devk510 * ztc ) 
     743         aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     744 
     745         ! CONVERT FROM DIFFERENT PH SCALES 
     746         total2free  = 1.0/(1.0 + zst/aks3s(ji)) 
     747         free2SWS    = 1. + zst/aks3s(ji) + zft/akf3s(ji) 
     748         total2SWS   = total2free * free2SWS 
     749         SWS2total   = 1.0 / total2SWS 
     750 
     751         ! Convert to total scale 
     752         ak1s(ji)  = ak1s(ji)  * SWS2total 
     753         ak2s(ji)  = ak2s(ji)  * SWS2total 
     754         akbs(ji)  = akbs(ji)  * SWS2total 
     755         akws(ji)  = akws(ji)  * SWS2total 
     756         ak1ps(ji) = ak1ps(ji) * SWS2total 
     757         ak2ps(ji) = ak2ps(ji) * SWS2total 
     758         ak3ps(ji) = ak3ps(ji) * SWS2total 
     759         aksis(ji) = aksis(ji) * SWS2total 
     760         akf3s(ji) = akf3s(ji) / total2free 
     761 
     762         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
     763         !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO 
     764         !        (P. 1285) AND BERNER (1976) 
     765         zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 
     766         zbuf2  = 0.5 * ( devk46 + devk56 * ztc ) 
     767         aksps(ji) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     768 
     769         ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 
     770         borats(ji)   = 0.0002414 * zcl / 10.811 
     771         sulfats(ji)  = zst 
     772         fluorids(ji) = zft 
     773 
     774         ! Iron and SIO3 saturation concentration from ... 
     775         sieqs(ji) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6 
     776      END DO 
     777      ! 
     778      IF( ln_timing )  CALL timing_stop('sed_chem_cst') 
     779      ! 
     780   END SUBROUTINE sed_chem_cst 
     781 
    572782 
    573783END MODULE sedchem 
Note: See TracChangeset for help on using the changeset viewer.