Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zche.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zche.F90	(revision 13295)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zche.F90	(revision 14086)
@@ -188,259 +188,238 @@
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
-         DO ji = 1, jpi
-            !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tempis(ji,jj,1) + 273.15
-            zt    = ztkel * 0.01
-            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
-            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
-            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
-            &       + 0.0047036e-4*ztkel**2)
-            chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
-            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
-            chemc(ji,jj,3) = 57.7 - 0.118*ztkel
-            !
-         END DO
-      END DO
+      DO_2D( 1, 1, 1, 1 )
+         !                             ! SET ABSOLUTE TEMPERATURE
+         ztkel = tempis(ji,jj,1) + 273.15
+         zt    = ztkel * 0.01
+         zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
+         !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
+         !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
+         zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
+         &       + 0.0047036e-4*ztkel**2)
+         chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
+         chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
+         chemc(ji,jj,3) = 57.7 - 0.118*ztkel
+      END_2D
 
       ! OXYGEN SOLUBILITY - DEEP OCEAN
       ! -------------------------------
-!CDIR NOVERRCHK
-      DO jk = 1, jpk
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-              ztkel = tempis(ji,jj,jk) + 273.15
-              zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
-              zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
-              ztgg2 = ztgg  * ztgg
-              ztgg3 = ztgg2 * ztgg
-              ztgg4 = ztgg3 * ztgg
-              ztgg5 = ztgg4 * ztgg
-
-              zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
-              &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
-              &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
-              &       - 3.11680e-7 * zsal2
-              chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
-            END DO
-          END DO
-        END DO
+      DO_3D( 1, 1, 1, 1, 1, jpk )
+         ztkel = tempis(ji,jj,jk) + 273.15
+         zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
+         zsal2 = zsal * zsal
+         ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+         ztgg2 = ztgg  * ztgg
+         ztgg3 = ztgg2 * ztgg
+         ztgg4 = ztgg3 * ztgg
+         ztgg5 = ztgg4 * ztgg
+
+         zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
+         &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
+         &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
+         &       - 3.11680e-7 * zsal2
+         chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
+      END_3D
 
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
-!CDIR NOVERRCHK
-      DO jk = 1, jpk
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-
-               ! SET PRESSION ACCORDING TO SAUNDER (1980)
-               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
-               zc1 = 5.92E-3 + zplat**2 * 5.25E-3
-               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6
-               zpres = zpres / 10.0
-
-               ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tempis(ji,jj,jk) + 273.15
-               zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
-               zsqrt  = SQRT( zsal )
-               zsal15  = zsqrt * zsal
-               zlogt  = LOG( ztkel )
-               ztr    = 1. / ztkel
-               zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
-               zis2   = zis * zis
-               zisqrt = SQRT( zis )
-               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
-
-               ! CHLORINITY (WOOSTER ET AL., 1969)
-               zcl     = zsal / 1.80655
-
-               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
-               zst     = 0.14 * zcl /96.062
-
-               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
-               zft     = 0.000067 * zcl /18.9984
-
-               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
-               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
-               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
-               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
-               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
-               &         + LOG(1.0 - 0.001005 * zsal))
-
-               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
-               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
-               &         + LOG(1.0d0 - 0.001005d0*zsal)            &
-               &         + LOG(1.0d0 + zst/zcks))
-
-               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
-               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
-               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
-               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
-               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
-               &      * zlogt + 0.053105*zsqrt*ztkel
-
-               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
-               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
-               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
-                  - 0.011555*zsal + 0.0001152*zsal*zsal)
-               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
-                  - 0.01781*zsal + 0.0001122*zsal*zsal)
-
-               ! PKW (H2O) (MILLERO, 1995) from composite data
-               zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
-                         - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
-
-               ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
-              zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
-              &          + (-106.736*ztr + 0.69171) * zsqrt       &
-              &          + (-0.65643*ztr - 0.01844) * zsal
-
-              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
-              &          + (-160.340*ztr + 1.3566)*zsqrt          &
-              &          + (0.37335*ztr - 0.05778)*zsal
-
-              zck3p    = -3070.75*ztr - 18.126                    &
-              &          + (17.27039*ztr + 2.81197) * zsqrt       &
-              &          + (-44.99486*ztr - 0.09984) * zsal 
-
-              ! CONSTANT FOR SILICATE, MILLERO (1995)
-              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
-              &          + (-458.79*ztr + 3.5913) * zisqrt       &
-              &          + (188.74*ztr - 1.5998) * zis           &
-              &          + (-12.1652*ztr + 0.07871) * zis2       &
-              &          + LOG(1.0 - 0.001005*zsal)
-
-               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
-               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
-               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
-                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
-                  &      - 0.07711*zsal + 0.0041249*zsal15
-
-               ! CONVERT FROM DIFFERENT PH SCALES
-               total2free  = 1.0/(1.0 + zst/zcks)
-               free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
-               total2SWS   = total2free * free2SWS
-               SWS2total   = 1.0 / total2SWS
-
-               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
-               zak1    = 10**(zck1) * total2SWS
-               zak2    = 10**(zck2) * total2SWS
-               zakb    = EXP( zckb ) * total2SWS
-               zakw    = EXP( zckw )
-               zaksp1  = 10**(zaksp0)
-               zak1p   = exp( zck1p )
-               zak2p   = exp( zck2p )
-               zak3p   = exp( zck3p )
-               zaksi   = exp( zcksi )
-               zckf    = zckf * total2SWS
-
-               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
-               !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
-               !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
-               !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
-               !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
-               !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
-               !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
-               !        & GIESKES (1970), P. 1285-1286 (THE SMALL
-               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
-               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
-               zcpexp  = zpres / (rgas*ztkel)
-               zcpexp2 = zpres * zcpexp
-
-               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
-               !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
-               !        (CF. BROECKER ET AL., 1982)
-
-               zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
-               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
-               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
-               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
-               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
-               aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
-               akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
-               ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
-               ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
-               ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
-               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               ! CONVERT FROM DIFFERENT PH SCALES
-               total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
-               free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
-               total2SWS   = total2free * free2SWS
-               SWS2total   = 1.0 / total2SWS
-
-               ! Convert to total scale
-               ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
-               ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
-               akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
-               akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
-               ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
-               ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
-               ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
-               aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
-               akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
-
-               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
-               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
-               !        (P. 1285) AND BERNER (1976)
-               zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
-               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
-               borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
-               sulfat(ji,jj,jk) = zst
-               fluorid(ji,jj,jk) = zft 
-
-               ! Iron and SIO3 saturation concentration from ...
-               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
-               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel )
-
-               ! Liu and Millero (1999) only valid 5 - 50 degC
-               ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
-               fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1)
-               fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 )
-               fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
-               fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
-               fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 )
-            END DO
-         END DO
-      END DO
+      DO_3D( 1, 1, 1, 1, 1, jpk )
+          ! SET PRESSION ACCORDING TO SAUNDER (1980)
+          zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
+          zc1 = 5.92E-3 + zplat**2 * 5.25E-3
+          zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6
+          zpres = zpres / 10.0
+
+          ! SET ABSOLUTE TEMPERATURE
+          ztkel   = tempis(ji,jj,jk) + 273.15
+          zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+          zsqrt  = SQRT( zsal )
+          zsal15  = zsqrt * zsal
+          zlogt  = LOG( ztkel )
+          ztr    = 1. / ztkel
+          zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
+          zis2   = zis * zis
+          zisqrt = SQRT( zis )
+          ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+
+          ! CHLORINITY (WOOSTER ET AL., 1969)
+          zcl     = zsal / 1.80655
+
+          ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
+          zst     = 0.14 * zcl /96.062
+
+          ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
+          zft     = 0.000067 * zcl /18.9984
+
+          ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
+          zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
+          &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
+          &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
+          &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
+          &         + LOG(1.0 - 0.001005 * zsal))
+
+          ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
+          zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
+          &         + LOG(1.0d0 - 0.001005d0*zsal)            &
+          &         + LOG(1.0d0 + zst/zcks))
+
+          ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
+          zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
+          &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
+          &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
+          &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
+          &      * zlogt + 0.053105*zsqrt*ztkel
+
+          ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
+          ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+          zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+             - 0.011555*zsal + 0.0001152*zsal*zsal)
+          zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+             - 0.01781*zsal + 0.0001122*zsal*zsal)
+
+          ! PKW (H2O) (MILLERO, 1995) from composite data
+          zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
+                    - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
+
+          ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
+         zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
+         &          + (-106.736*ztr + 0.69171) * zsqrt       &
+         &          + (-0.65643*ztr - 0.01844) * zsal
+
+         zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
+         &          + (-160.340*ztr + 1.3566)*zsqrt          &
+         &          + (0.37335*ztr - 0.05778)*zsal
+
+         zck3p    = -3070.75*ztr - 18.126                    &
+         &          + (17.27039*ztr + 2.81197) * zsqrt       &
+         &          + (-44.99486*ztr - 0.09984) * zsal 
+
+         ! CONSTANT FOR SILICATE, MILLERO (1995)
+         zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
+         &          + (-458.79*ztr + 3.5913) * zisqrt       &
+         &          + (188.74*ztr - 1.5998) * zis           &
+         &          + (-12.1652*ztr + 0.07871) * zis2       &
+         &          + LOG(1.0 - 0.001005*zsal)
+
+          ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
+          !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
+          zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+             &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+             &      - 0.07711*zsal + 0.0041249*zsal15
+
+          ! CONVERT FROM DIFFERENT PH SCALES
+          total2free  = 1.0/(1.0 + zst/zcks)
+          free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
+          total2SWS   = total2free * free2SWS
+          SWS2total   = 1.0 / total2SWS
+
+          ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+          zak1    = 10**(zck1) * total2SWS
+          zak2    = 10**(zck2) * total2SWS
+          zakb    = EXP( zckb ) * total2SWS
+          zakw    = EXP( zckw )
+          zaksp1  = 10**(zaksp0)
+          zak1p   = exp( zck1p )
+          zak2p   = exp( zck2p )
+          zak3p   = exp( zck3p )
+          zaksi   = exp( zcksi )
+          zckf    = zckf * total2SWS
+
+          ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
+          !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
+          !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
+          !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
+          !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
+          !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
+          !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
+          !        & GIESKES (1970), P. 1285-1286 (THE SMALL
+          !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
+          !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
+          zcpexp  = zpres / (rgas*ztkel)
+          zcpexp2 = zpres * zcpexp
+
+          ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+          !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
+          !        (CF. BROECKER ET AL., 1982)
+
+          zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
+          ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
+          ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
+          akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
+          akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
+          aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+          akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
+          ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
+          ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
+          ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
+          aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          ! CONVERT FROM DIFFERENT PH SCALES
+          total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
+          free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
+          total2SWS   = total2free * free2SWS
+          SWS2total   = 1.0 / total2SWS
+
+          ! Convert to total scale
+          ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
+          ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
+          akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
+          akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
+          ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
+          ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
+          ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
+          aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
+          akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
+
+          ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
+          !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
+          !        (P. 1285) AND BERNER (1976)
+          zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
+          aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
+          borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
+          sulfat(ji,jj,jk) = zst
+          fluorid(ji,jj,jk) = zft 
+
+          ! Iron and SIO3 saturation concentration from ...
+          sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
+          fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 
+          ! Liu and Millero (1999) only valid 5 - 50 degC
+          ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
+          fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1)
+          fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 )
+          fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
+          fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
+          fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 )
+      END_3D
       !
       IF( ln_timing )  CALL timing_stop('p4z_che')
