Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F	(revision 772)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90	(revision 775)
@@ -1,310 +1,258 @@
-
-CCC $Header$ 
-CCC  TOP 1.0 , LOCEAN-IPSL (2005) 
-C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-C ---------------------------------------------------------------------------
-CDIR$ LIST
-      SUBROUTINE p4zche
-#if defined key_top && defined key_pisces
-CCC---------------------------------------------------------------------
-CCC
-CCC          ROUTINE p4zche : PISCES MODEL
-CCC          *****************************
-CCC
-CCC     PURPOSE.
-CCC     --------
-CCC          *P4ZCHE* : Sea water chemistry computed following OCMIP protocol
-CCC
-CCC
-CC     EXTERNALS.
-CC     ----------
-CC          rhop
-CC
-CC   MODIFICATIONS:
-CC   --------------
-CC      original :      1988 E. Maier-Reimer
-CC      additions :     1998 O. Aumont
-CC      modifications : 1999 C. Le Quere
-CC      modifications : 2004 O. Aumont
-CC      modifications : 2006 R. Gangsto
-CC----------------------------------------------------------------------
-CC parameters and commons
-CC ======================
-CDIR$ nolist
-      USE oce_trc
-      USE trp_trc
-      USE sms
-      IMPLICIT NONE
+MODULE p4zche
+   !!======================================================================
+   !!                         ***  MODULE p4zche  ***
+   !! TOP :   PISCES Sea water chemistry computed following OCMIP protocol
+   !!======================================================================
+   !! History :    -   !  1988     (E. Maier-Reimer)  Original code
+   !!              -   !  1998     (O. Aumont)  addition
+   !!              -   !  1999     (C. Le Quere)  modification
+   !!             1.0  !  2004     (O. Aumont)  modification
+   !!              -   !  2006     (R. Gangsto)  modification
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!----------------------------------------------------------------------
+#if defined key_pisces
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !!   p4z_che        :  Sea water chemistry computed following OCMIP protocol
+   !!----------------------------------------------------------------------
+   USE oce_trc         !
+   USE trp_trc         !
+   USE sms             ! 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_che    ! called in p4zprg.F90
+
+   !!* Substitution
 #include "domzgr_substitute.h90"
-CDIR$ list
-CC----------------------------------------------------------------------
-CC local declarations
-CC ==================
-C
-      INTEGER ji, jj, jk
-      REAL tkel, sal,  qtt, zbuf1, zbuf2
-      REAL pres, tc, cl, cpexp, cek0, oxy, cpexp2
-      REAL zsqrt, ztr, zlogt, cek1
-      REAL zqtt, qtt2, sal15, zis, zis2, zisqrt
-      REAL ckb, ck1, ck2, ckw, ak1, ak2, akb, aksp0, akw
-      REAL ckp1, ckp2, ckp3, cksi, akp1, akp2, akp3, aksi
-      REAL st, ft, cks, ckf, aks, akf, aksp1
-
-C
-C* 1. CHEMICAL CONSTANTS - SURFACE LAYER
-C ---------------------------------------
-C
-      DO jj = 1,jpj
-        DO ji = 1,jpi
-C
-C* 1.1 SET ABSOLUTE TEMPERATURE
-C ------------------------------
-C
-          tkel = tn(ji,jj,1)+273.16
-          qtt = tkel*0.01
-          qtt2=qtt*qtt
-          sal = sn(ji,jj,1) + (1.-tmask(ji,jj,1))*35.
-          zqtt=log(qtt)
-C
-C* 1.2 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
-C      AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-C -------------------------------------------------------
-C
-          cek0 = c00+c01/qtt+c02*zqtt+sal*(c03+c04*qtt+c05*qtt2)
-          cek1 = ca0+ca1/qtt+ca2*zqtt+ca3*qtt2+sal*(ca4
-     &      +ca5*qtt+ca6*qtt2)
-C
-C* 1.3 LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970)
-C ------------------------------------------------------------
-C
-          oxy = ox0+ox1/qtt+ox2*zqtt+sal*(ox3+ox4*qtt+ox5*qtt2)
-C
-C* 1.4 SET SOLUBILITIES OF O2 AND CO2
-C -----------------------------------
-C
-          chemc(ji,jj,1) = exp(cek0)*1.E-6*rhop(ji,jj,1)/1000.
-          chemc(ji,jj,2) = exp(oxy)*oxyco
-          chemc(ji,jj,3) = exp(cek1)*1.E-6*rhop(ji,jj,1)/1000.
-C
-        ENDDO
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Header:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE p4z_che
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_che  ***
+      !!
+      !! ** Purpose :   Sea water chemistry computed following OCMIP protocol
+      !!
+      !! ** Method  : - ...
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   ztkel, zsal , zqtt  , zbuf1 , zbuf2
+      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zcek0, zoxy  , zcpexp2
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
+      REAL(wp) ::   zlqtt, zqtt2, zsal15, zis   , zis2 , zisqrt
+      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
+      REAL(wp) ::   zckp1, zckp2, zckp3 , zcksi , zakp1, zakp2 , zakp3, zaksi
+      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaks , zakf  , zaksp1
+      !!---------------------------------------------------------------------
+
+      ! CHEMICAL CONSTANTS - SURFACE LAYER
+      ! ----------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            !                             ! SET ABSOLUTE TEMPERATURE
+            ztkel = tn(ji,jj,1) + 273.16
+            zqtt  = ztkel * 0.01
+            zqtt2 = zqtt * zqtt
+            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
+            zlqtt = LOG( zqtt )
+
+            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
+            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
+            zcek0 = c00 + c01 / zqtt + c02 * zlqtt + zsal * ( c03 + c04 * zqtt + c05 * zqtt2 )
+            zcek1 = ca0 + ca1 / zqtt + ca2 * zlqtt + ca3 * zqtt2 + zsal*( ca4 + ca5 * zqtt + ca6 * zqtt2 )
+
+            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970)
+            zoxy  = ox0 + ox1 / zqtt + ox2 * zlqtt + zsal * ( ox3 + ox4 * zqtt + ox5 * zqtt2 )
+
+            !                             ! SET SOLUBILITIES OF O2 AND CO2
+            chemc(ji,jj,1) = EXP( zcek0 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
+            chemc(ji,jj,2) = EXP( zoxy  ) * oxyco
+            chemc(ji,jj,3) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
+
+         END DO
       END DO
-C
-C* 2 CHEMICAL CONSTANTS - DEEP OCEAN
-C -------------------------------------
-C
-      DO jk = 1,jpk
-        DO jj = 1,jpj
-          DO ji = 1,jpi
-C
-C* 2.1 SET PRESSION
-C -----------------
-C
-            pres = 1.025e-1*fsdept(ji,jj,jk)
-C
-C* 2.2 SET ABSOLUTE TEMPERATURE
-C ------------------------------
-C
-            tkel   = tn(ji,jj,jk)+273.16
-            qtt    = tkel*0.01
-            sal    = sn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*35.
-            zsqrt  = sqrt(sal)
-            sal15  = zsqrt*sal
-            zlogt  = log(tkel)
-            ztr    = 1./tkel
-            zis    = 19.924*sal/(1000.-1.005*sal)
-            zis2   = zis*zis
-            zisqrt = sqrt(zis)
-            tc = tn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*20.
-C
-C* 2.3 CHLORINITY (WOOSTER ET AL., 1969)
-C ---------------------------------------
-C
-            cl = sal*salchl
-C
-C* 2.4 TOTAL SULFATE CONCENTR. [MOLES/kg soln]
-C --------------------------------------------
-C
-            st = st1*cl*st2
-C
-C* 2.5 TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
-C ---------------------------------------------
-C
-            ft = ft1*cl*ft2
-C
-C* 2.6 DISSOCIATION CONSTANT FOR SULFATES
-C on free H scale (Dickson 1990)
-C -------------------------------------------------------
-C
-            cks=exp(ks1*ztr+ks0+ks2*zlogt+(ks3*ztr+ks4+ks5*zlogt)
-     &      *zisqrt+(ks6*ztr+ks7+ks8*zlogt)*zis+ks9*ztr*zis
-     &      *zisqrt+ks10*ztr*zis2+log(ks11+ks12*sal))
-C
-C* 2.7 DISSOCIATION CONSTANT FOR FLUORIDES
-C on free H scale (Dickson and Riley 79)
-C -------------------------------------------------------
-C
-            ckf=exp(kf1*ztr+kf0+kf2*zisqrt+log(kf3+kf4*sal))
-
-C
-C* 2.4 DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
-C -------------------------------------------------------
-C
-            ckb = (cb0+cb1*zsqrt+cb2*sal+cb3*sal15+cb4*sal*sal)*ztr
-     &          +(cb5+cb6*zsqrt+cb7*sal)+
-     &          (cb8+cb9*zsqrt+cb10*sal)*zlogt+cb11*zsqrt*tkel
-     &          +log((1.+st/cks+ft/ckf)/(1.+st/cks))
-            ck1 = c10*ztr+c11+c12*zlogt+c13*sal+c14*sal**2
-            ck2 = c20*ztr+c21+c22*sal+c23*sal**2
-C
-C* 2.5 PKW (H2O) (DICKSON AND RILEY, 1979)
-C -----------------------------------------
-C
-            ckw = cw0*ztr+cw1+cw2*zlogt+(cw3*ztr+cw4+cw5*zlogt)*
-     &          zsqrt+cw6*sal
-
-C
-C
-C* 2.10 DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale)
-C  ---------------------------------------------------------------------
-C
-            ckp1 = cp10+cp11*ztr+cp12*zlogt+zsqrt*(cp13*ztr
-     &           +cp14)+sal*(cp15*ztr+cp16)
-            ckp2 = cp20+cp21*ztr+cp22*zlogt+zsqrt*(cp23*ztr
-     &           +cp24)+sal*(cp25*ztr+cp26)
-            ckp3 = cp30+cp31*ztr+zsqrt*(cp32*ztr
-     &           +cp33)+sal*(cp34*ztr+cp35)
-            cksi = cs10+cs11*ztr+cs12*zlogt+zisqrt*(cs13*ztr
-     &           +cs14)+zis*(cs15*ztr+cs16)+zis2*(cs17*ztr
-     &           +cs18)+log(1.+cs19*sal)
-     &           +log(cs20+cs21*sal)
-
-C
-C*2.7 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
-C       (S=27-43, T=2-25 DEG C) AT pres =0 (ATMOSPH. PRESSURE)
-C       (MUCCI 1983)
-C -------------------------------------------------------------
-C
-            aksp0 = akcc1+akcc2*tkel+akcc3*ztr+akcc4*log10(tkel)+
-     &      (akcc5+akcc6*tkel+
-     &      akcc7*ztr)*zsqrt+akcc8*sal+akcc9*sal15
-
-
-C
-C* 2.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
-C -----------------------------------------------------------------
-C
-            ak1   = 10**(ck1)
-            ak2   = 10**(ck2)
-            akb   = exp(ckb)
-            akp1  = exp(ckp1)
-            akp2  = exp(ckp2)
-            akp3  = exp(ckp3)
-            aksi  = exp(cksi)
-            akw   = exp(ckw)
-            aksp1 = 10**(aksp0)
-            aks   = exp(cks)
-            akf   = exp(ckf)
-
-
-C
-C* 2.8 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
-C        (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE
-C        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
-C        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
-C        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
-C        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
-C        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
-C        AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL
-C        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
-C        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON
-C        P. 1285))
-C -----------------------------------------------------------
-C
-            cpexp = pres /(rgas*tkel)
-            cpexp2 = pres * pres/(rgas*tkel)
-C
-C* 2.9 KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
-C        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
-C        (CF. BROECKER ET AL., 1982)
-C --------------------------------------------------------
-C
-            zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc)
-            zbuf2 = 0.5*(devk4(3)+devk5(3)*tc)
-            akb3(ji,jj,jk) = akb*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(1)+devk2(1)*tc+devk3(1)*tc*tc)
-            zbuf2 = 0.5*(devk4(1)+devk5(1)*tc)
-            ak13(ji,jj,jk) = ak1*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(2)+devk2(2)*tc+devk3(2)*tc*tc)
-            zbuf2 = 0.5*(devk4(2)+devk5(2)*tc)
-            ak23(ji,jj,jk) = ak2*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(4)+devk2(4)*tc+devk3(4)*tc*tc)
-            zbuf2 = 0.5*(devk4(4)+devk5(4)*tc)
-            akp13(ji,jj,jk) = akp1*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(5)+devk2(5)*tc+devk3(5)*tc*tc)
-            zbuf2 = 0.5*(devk4(5)+devk5(5)*tc)
-            akp23(ji,jj,jk) = akp2*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(6)+devk2(6)*tc+devk3(6)*tc*tc)
-            zbuf2 = 0.5*(devk4(6)+devk5(6)*tc)
-            akp33(ji,jj,jk) = akp3*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(7)+devk2(7)*tc+devk3(7)*tc*tc)
-            zbuf2 = 0.5*(devk4(7)+devk5(7)*tc)
-            akw3(ji,jj,jk) = akw*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-C  Ksi
-C            aksi3(ji,jj,jk) = aksi
-C
-C  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys)
-C  "deltaVsi and deltaKsi have been estimated from the value of boric acid"
-C
-            zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc)
-            zbuf2 = 0.5*(devk4(3)+devk5(3)*tc)
-            aksi3(ji,jj,jk) = aksi*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-C
-C
-C* 2.15 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
-C        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
-C        (P. 1285) AND BERNER (1976)
-C -------------------------------------------------
-
-            zbuf1 = -(devk1(8)+devk2(8)*tc+devk3(8)*tc*tc)
-            zbuf2 = 0.5*(devk4(8)+devk5(8)*tc)
-            aksp(ji,jj,jk) = aksp1*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-C  Pressure correction for sulfate and fluoride
-C
-            zbuf1 = -(devk1(9)+devk2(9)*tc+devk3(9)*tc*tc)
-            zbuf2 = 0.5*(devk4(9)+devk5(9)*tc)
-            aks3(ji,jj,jk) = aks*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-            zbuf1 = -(devk1(10)+devk2(10)*tc+devk3(10)*tc*tc)
-            zbuf2 = 0.5*(devk4(10)+devk5(10)*tc)
-            akf3(ji,jj,jk) = akf*exp(zbuf1*cpexp+zbuf2*cpexp2)
-
-
-C
-C* 2.11 TOTAL BORATE CONCENTR. [MOLES/L]
-C --------------------------------------
-C
-            borat(ji,jj,jk) = bor1*cl*bor2
-C
-C  2.12 Iron and SIO3 saturation concentration from ...
-C  ----------------------------------------------------
-C
-         sio3eq(ji,jj,jk)=exp(log(10.)*(6.44-968./tkel))*1E-6
-         fekeq(ji,jj,jk)=10**(17.27-1565.7/(273.15+tc))
-C
-          ENDDO
-        ENDDO
+
+      ! CHEMICAL CONSTANTS - DEEP OCEAN
+      ! -------------------------------
+
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               ! SET PRESSION
+               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+
+               ! SET ABSOLUTE TEMPERATURE
+               ztkel   = tn(ji,jj,jk) + 273.16
+               zqtt    = ztkel * 0.01
+               zsal    = sn(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     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+
+               ! CHLORINITY (WOOSTER ET AL., 1969)
+               zcl     = zsal * salchl
+
+               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
+               zst     = st1 * zcl * st2
+
+               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
+               zft     = ft1 * zcl * ft2
+
+               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
+               zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           &
+                  &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   &
+                  &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      &
+                  &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  )
+
+               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
+               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  )
+
+               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
+               zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   &
+                  &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            &
+                  &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             &
+                  &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  )
+!!gm zsal**2 to be replaced by a *...
+               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal**2
+               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
+
+               ! PKW (H2O) (DICKSON AND RILEY, 1979)
+               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
+
+               ! DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale)
+               zckp1 = cp10 + cp11 * ztr + cp12 * zlogt + zsqrt * ( cp13 * ztr + cp14 ) + zsal * ( cp15 * ztr + cp16 )
+               zckp2 = cp20 + cp21 * ztr + cp22 * zlogt + zsqrt * ( cp23 * ztr + cp24 ) + zsal * ( cp25 * ztr + cp26 )
+               zckp3 = cp30 + cp31 * ztr                + zsqrt * ( cp32 * ztr + cp33 ) + zsal * ( cp34 * ztr + cp35 )
+               zcksi = cs10 + cs11 * ztr + cs12 * zlogt + zisqrt* ( cs13 * ztr + cs14 ) + zis * ( cs15 * ztr + cs16 )  &
+                  &                                    + zis2  * ( cs17 * ztr + cs18 ) + LOG( 1.   + cs19 * zsal )     &
+                  &                                                                    + LOG( cs20 + cs21 * 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  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
+                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
+
+               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+               zak1    = 10**(zck1)
+               zak2    = 10**(zck2)
+               zakb    = EXP( zckb  )
+               zakp1   = EXP( zckp1 )
+               zakp2   = EXP( zckp2 )
+               zakp3   = EXP( zckp3 )
+               zaksi   = EXP( zcksi )
+               zakw    = EXP( zckw )
+               zaksp1  = 10**(zaksp0)
+               zaks    = exp( zcks )
+               zakf    = exp( zckf )
+
+               ! 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 * zpres/(rgas*ztkel)
+
+               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+               !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
+               !        (CF. BROECKER ET AL., 1982)
+               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
+               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  = -(devk1(1)+devk2(1)*ztc+devk3(1)*ztc*ztc)
+               zbuf2  = 0.5*(devk4(1)+devk5(1)*ztc)
+               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(2) + devk2(2) * ztc + devk3(2) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(2) + devk5(2) * ztc )
+               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(4) + devk2(4) * ztc + devk3(4) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(4) + devk5(4) * ztc )
+               akp13(ji,jj,jk) = zakp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(5) + devk2(5) * ztc + devk3(5) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(5) + devk5(5) * ztc )
+               akp23(ji,jj,jk) = zakp2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(6) + devk2(6) * ztc + devk3(6) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(6) + devk5(6) * ztc )
+               akp33(ji,jj,jk) = zakp3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(7) + devk2(7) * ztc + devk3(7) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(7) + devk5(7) * ztc )
+               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               !  Ksi
+               !            aksi3(ji,jj,jk) = zaksi
+               !
+               !  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys)
+               !  "deltaVsi and deltaKsi have been estimated from the value of boric acid"
+               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
+               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
+               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
+               !        (P. 1285) AND BERNER (1976)
+               zbuf1  =     - ( devk1(8) + devk2(8) * ztc + devk3(8) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(8) + devk5(8) * ztc )
+               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               !  Pressure correction for sulfate and fluoride
+               zbuf1  =     - ( devk1(9) + devk2(9) * ztc + devk3(9) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(9) + devk5(9) * ztc )
+               aks3(ji,jj,jk) = zaks   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk1(10) + devk2(10) * ztc + devk3(10) * ztc * ztc )
+               zbuf2  = 0.5 * ( devk4(10) + devk5(10) * ztc )
+               akf3(ji,jj,jk) = zakf   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               ! TOTAL BORATE CONCENTR. [MOLES/L]
+               borat(ji,jj,jk) = bor1 * zcl * bor2
+
+               ! 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 / ( 273.15 + ztc ) )
+
+            END DO
+         END DO
       END DO
-C
-#endif
-C
-      RETURN
-      END
+      !
+   END SUBROUTINE p4z_che
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No PISCES bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE p4z_che( kt )                   ! Empty routine
+      INTEGER, INTENT( in ) ::   kt
+      WRITE(*,*) 'p4z_che: You should not have seen this print! error?', kt
+   END SUBROUTINE p4z_che
+#endif 
+
+   !!======================================================================
+END MODULE  p4zche
