Index: branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 9816)
+++ branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 9817)
@@ -32,5 +32,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
 
    REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm
@@ -39,15 +40,5 @@
    REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )  
 
-   REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium
-   REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983
-   REAL(wp) ::   akcc3  = 2839.319        
-   REAL(wp) ::   akcc4  =   71.595        
-   REAL(wp) ::   akcc5  =   -0.77712      
-   REAL(wp) ::   akcc6  =    0.00284263   
-   REAL(wp) ::   akcc7  =  178.34        
-   REAL(wp) ::   akcc8  =   -0.07711     
-   REAL(wp) ::   akcc9  =    0.0041249   
-
-   REAL(wp) ::   rgas   = 83.143         ! universal gas constants
+   REAL(wp) ::   rgas   = 83.14472       ! universal gas constants
    REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles
 
@@ -55,67 +46,9 @@
    REAL(wp) ::   bor2   = 1. / 10.82
 
-   REAL(wp) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm)
-   REAL(wp) ::   ca1    =  218.2968
-   REAL(wp) ::   ca2    =   90.9241
-   REAL(wp) ::   ca3    =   -1.47696
-   REAL(wp) ::   ca4    =    0.025695
-   REAL(wp) ::   ca5    =   -0.025225
-   REAL(wp) ::   ca6    =    0.0049867
-
-   REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)   
-   REAL(wp) ::   c11    =    62.008     
-   REAL(wp) ::   c12    =    -9.7944    
-   REAL(wp) ::   c13    =     0.0118     
-   REAL(wp) ::   c14    =    -0.000116
-
-   REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
-   REAL(wp) ::   c21    =    -4.777   
-   REAL(wp) ::   c22    =     0.0184   
-   REAL(wp) ::   c23    =    -0.000118
-
    REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate
    REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966)
-   REAL(wp) ::   ks0    =    141.328 
-   REAL(wp) ::   ks1    =  -4276.1  
-   REAL(wp) ::   ks2    =    -23.093
-   REAL(wp) ::   ks3    = -13856.  
-   REAL(wp) ::   ks4    =   324.57 
-   REAL(wp) ::   ks5    =   -47.986
-   REAL(wp) ::   ks6    =  35474. 
-   REAL(wp) ::   ks7    =   -771.54
-   REAL(wp) ::   ks8    =    114.723
-   REAL(wp) ::   ks9    =  -2698.  
-   REAL(wp) ::   ks10   =   1776. 
-   REAL(wp) ::   ks11   =      1.
-   REAL(wp) ::   ks12   =     -0.001005 
 
    REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides
    REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 )
-   REAL(wp) ::   kf0    =  -12.641    
-   REAL(wp) ::   kf1    = 1590.2    
-   REAL(wp) ::   kf2    =    1.525    
-   REAL(wp) ::   kf3    =    1.0     
-   REAL(wp) ::   kf4    =   -0.001005
-
-   REAL(wp) ::   cb0    = -8966.90      ! Coeff. for 1. dissoc. of boric acid 
-   REAL(wp) ::   cb1    = -2890.53      ! (Dickson and Goyet, 1994)
-   REAL(wp) ::   cb2    =   -77.942
-   REAL(wp) ::   cb3    =     1.728
-   REAL(wp) ::   cb4    =    -0.0996
-   REAL(wp) ::   cb5    =   148.0248
-   REAL(wp) ::   cb6    =   137.1942
-   REAL(wp) ::   cb7    =     1.62142
-   REAL(wp) ::   cb8    =   -24.4344
-   REAL(wp) ::   cb9    =   -25.085
-   REAL(wp) ::   cb10   =    -0.2474 
-   REAL(wp) ::   cb11   =     0.053105
-
-   REAL(wp) ::   cw0    = -13847.26     ! Coeff. for dissoc. of water (Dickson and Riley, 1979 )
-   REAL(wp) ::   cw1    =    148.9652  
-   REAL(wp) ::   cw2    =    -23.6521
-   REAL(wp) ::   cw3    =    118.67 
-   REAL(wp) ::   cw4    =     -5.977 
-   REAL(wp) ::   cw5    =      1.0495  
-   REAL(wp) ::   cw6    =     -0.01615
 
    !                                    ! volumetric solubility constants for o2 in ml/L  
@@ -185,6 +118,6 @@
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
-      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
-      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
       REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
@@ -193,4 +126,20 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_che')
       !
+      ! Computations of chemical constants require in situ temperature
+      ! Here a quite simple formulation is used to convert 
+      ! potential temperature to in situ temperature. The errors is less than 
+      ! 0.04°C relative to an exact computation
+      ! ---------------------------------------------------------------------
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zpres = fsdept(ji,jj,jk) / 1000.
+               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) )
+               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 )
+               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2
+            END DO
+         END DO
+      END DO
+      !
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
@@ -200,5 +149,5 @@
          DO ji = 1, jpi
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.16
+            ztkel = tempis(ji,jj,1) + 273.15
             zt    = ztkel * 0.01
             zt2   = zt * zt
@@ -208,17 +157,10 @@
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 )
-            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992)
-            ztgg  = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
-            ztgg2 = ztgg  * ztgg
-            ztgg3 = ztgg2 * ztgg
-            ztgg4 = ztgg3 * ztgg
-            ztgg5 = ztgg4 * ztgg
-            zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   &
-                   + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2
-
+            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
+            &       + 0.0047036e-4*ztkel**2)
             !                             ! SET SOLUBILITIES OF O2 AND CO2 
-            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm)
-            chemc(ji,jj,2) = ( EXP( zoxy  ) * o2atm ) * oxyco              ! mol/(L atm)
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm)
+            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
@@ -233,8 +175,8 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-              ztkel = tsn(ji,jj,jk,jp_tem) + 273.16
+              ztkel = tempis(ji,jj,jk) + 273.15
               zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
               ztgg2 = ztgg  * ztgg
               ztgg3 = ztgg2 * ztgg
@@ -259,9 +201,12 @@
             DO ji = 1, jpi
 
-               ! SET PRESSION
-               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+               ! 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*fsdept(ji,jj,jk)))) / 4.42E-6
+               zpres = zpres / 10.0
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16
+               ztkel   = tempis(ji,jj,jk) + 273.15
                zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
@@ -272,5 +217,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
@@ -284,29 +229,42 @@
 
                ! 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 )  )
+               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))
+               !
+               aphscale(ji,jj,jk) = ( 1. + zst / zcks )
 
                ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
-               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  )
+               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    = ( 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 )  )
-
-               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal
-               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
+               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) (DICKSON AND RILEY, 1979)
-               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
-
+               zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    & 
+               &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        &
+               &     * zsqrt - 0.01615 * 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
+               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
 
                ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
@@ -378,5 +336,6 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc )
+      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   &
+      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc )
       !
       IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.')
@@ -396,3 +355,3 @@
 
    !!======================================================================
-END MODULE  p4zche
+END MODULE p4zche
