/[MITgcm]/MITgcm_contrib/darwin2/pkg/darwin/dic_carbon_chem.F
ViewVC logotype

Diff of /MITgcm_contrib/darwin2/pkg/darwin/dic_carbon_chem.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by stephd, Mon May 9 20:17:39 2011 UTC revision 1.4 by stephd, Wed Oct 9 15:37:05 2013 UTC
# Line 634  c======================================= Line 634  c=======================================
634  CStartOfInterFace  CStartOfInterFace
635        SUBROUTINE CARBON_COEFFS(        SUBROUTINE CARBON_COEFFS(
636       I                   ttemp,stemp,       I                   ttemp,stemp,
637       I                   bi,bj,iMin,iMax,jMin,jMax,myThid)       I                   bi,bj,iMin,iMax,jMin,jMax,
638         I                   kLevel, myThid)
639  C  C
640  C     /==========================================================\  C     /==========================================================\
641  C     | SUBROUTINE CARBON_COEFFS                                 |  C     | SUBROUTINE CARBON_COEFFS                                 |
# Line 642  C     | determine coefficients for surfa Line 643  C     | determine coefficients for surfa
643  C     | adapted from OCMIP2:  SUBROUTINE CO2CALC                 |  C     | adapted from OCMIP2:  SUBROUTINE CO2CALC                 |
644  C     | mick follows, oct 1999                                   |  C     | mick follows, oct 1999                                   |
645  c     | minor changes to tidy, swd aug 2002                      |  c     | minor changes to tidy, swd aug 2002                      |
646    c     | MODIFIED FOR PRESSURE DEPENDENCE                         |
647    c     | Karsten Friis and Mick Follows 2004                      |
648    c     | added 2013 (steph)                                       |
649  C     \==========================================================/  C     \==========================================================/
650  C INPUT  C INPUT
651  C       diclocal = total inorganic carbon (mol/m^3)  C       diclocal = total inorganic carbon (mol/m^3)
# Line 663  c       given in mol/m^3; output argumen Line 667  c       given in mol/m^3; output argumen
667  c       are likewise be in mol/m^3.  c       are likewise be in mol/m^3.
668  C  C
669  C Apr 2011: fix vapour bug (following Bennington)  C Apr 2011: fix vapour bug (following Bennington)
670    C Oct 2013: c NOW INCLUDES:
671    c PRESSURE DEPENDENCE of K1, K2, solubility product of calcite
672    c based on Takahashi, GEOSECS Atlantic Report, Vol. 1 (1981)
673  C--------------------------------------------------------------------------  C--------------------------------------------------------------------------
674          IMPLICIT NONE          IMPLICIT NONE
675  C     == GLobal variables ==  C     == GLobal variables ==
# Line 681  C general Line 688  C general
688          _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)          _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
689          _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)          _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
690          INTEGER bi,bj,iMin,iMax,jMin,jMax          INTEGER bi,bj,iMin,iMax,jMin,jMax
691            INTEGER kLevel
692          INTEGER myThid          INTEGER myThid
693  CEndOfInterface  CEndOfInterface
694    
# Line 706  C LOCAL VARIABLES Line 714  C LOCAL VARIABLES
714          _RL  invtk          _RL  invtk
715          _RL  is          _RL  is
716          _RL  is2          _RL  is2
717    c add pressure dependency
718            _RL  bdepth
719            _RL  cdepth
720            _RL  pressc
721    c calcite stuff
722            _RL  Ksp_T_Calc
723            _RL  xvalue
724            _RL  zdum
725            _RL  tmpa1
726            _RL  tmpa2
727            _RL  tmpa3
728            _RL  logKspc
729            _RL  dv
730            _RL  dk
731            _RL  pfactor
732            _RL  bigR
733  c add Bennington  c add Bennington
734          _RL  P1atm          _RL  P1atm
735          _RL  Rgas          _RL  Rgas
# Line 715  c add Bennington Line 739  c add Bennington
739          _RL  B          _RL  B
740          INTEGER i          INTEGER i
741          INTEGER j          INTEGER j
742            INTEGER k
743    
744  C.....................................................................  C.....................................................................
745  C OCMIP note:  C OCMIP note:
# Line 727  C "Handbook of Methods C for the Analysi Line 752  C "Handbook of Methods C for the Analysi
752  C the Carbon Dioxide System in Seawater", DOE, 1994 (SOP No. 3, p25-26).  C the Carbon Dioxide System in Seawater", DOE, 1994 (SOP No. 3, p25-26).
753  C....................................................................  C....................................................................
754    
755    c determine pressure (bar) from depth
756    c 1 BAR at z=0m (atmos pressure)
757    c use UPPER surface of cell so top layer pressure = 0 bar
758    c for surface exchange coeffs
759    
760    c if surface, calculate at interface pressure,
761    c else calculate at mid-depth pressure
762            if (Klevel.gt.1) then
763             bdepth = 0.0d0
764             cdepth = 0.0d0
765             pressc = 1.01325 _d 0
766             do k = 1,Klevel
767                 cdepth = bdepth + 0.5d0*drF(k)
768                 bdepth = bdepth + drF(k)
769                 pressc = 1.0d0 + 0.1d0*cdepth
770             end do
771            else
772              pressc = 1.01325 _d 0
773            endif
774    
775          do i=imin,imax          do i=imin,imax
776           do j=jmin,jmax           do j=jmin,jmax
777            if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then            if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then
# Line 750  C added by Val Bennington Nov 2010 Line 795  C added by Val Bennington Nov 2010
795  C Fugacity Factor needed for non-ideality in ocean  C Fugacity Factor needed for non-ideality in ocean
796  C ff used for atmospheric correction for water vapor and pressure  C ff used for atmospheric correction for water vapor and pressure
797  C Weiss (1974) Marine Chemistry  C Weiss (1974) Marine Chemistry
798             P1atm = 1.01325 _d 0 ! bars             P1atm = pressc ! bars
799             Rgas = 83.1451 _d 0 ! bar*cm3/(mol*K)             Rgas = 83.1451 _d 0 ! bar*cm3/(mol*K)
800             RT = Rgas*tk             RT = Rgas*tk
801             delta = (57.7 _d 0 - 0.118 _d 0*tk)             delta = (57.7 _d 0 - 0.118 _d 0*tk)
# Line 779  C Millero p.664 (1995) using Mehrbach et Line 824  C Millero p.664 (1995) using Mehrbach et
824       &          0.0118 _d 0 * s + 0.000116 _d 0*s2))       &          0.0118 _d 0 * s + 0.000116 _d 0*s2))
825             ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*invtk+ 4.777 _d 0-             ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*invtk+ 4.777 _d 0-
826       &          0.0184 _d 0*s + 0.000118 _d 0*s2))       &          0.0184 _d 0*s + 0.000118 _d 0*s2))
827    C
828    C NOW PRESSURE DEPENDENCE:
829    c Following Takahashi (1981) GEOSECS report - quoting Culberson and
830    c Pytkowicz (1968)
831               if (kLevel.gt.1) then
832    c pressc = pressure in bars
833               ak1(i,j,bi,bj) = ak1(i,j,bi,bj)*
834         &             exp( (24.2-0.085*t)*(pressc-1.0)/(83.143*tk) )
835    c FIRST GO FOR K2: According to GEOSECS (1982) report
836    c          ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*
837    c    &             exp( (26.4-0.040*t)*(pressc-1.0)/(83.143*tk) )
838    c SECOND GO FOR K2: corrected coeff according to CO2sys documentation
839    c                   E. Lewis and D. Wallace (1998) ORNL/CDIAC-105
840               ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*
841         &             exp( (16.4-0.040*t)*(pressc-1.0)/(83.143*tk) )
842               endif
843  C------------------------------------------------------------------------  C------------------------------------------------------------------------
844  C kb = [H][BO2]/[HBO2]  C kb = [H][BO2]/[HBO2]
845  C Millero p.669 (1995) using data from dickson (1990)  C Millero p.669 (1995) using data from dickson (1990)
# Line 787  C Millero p.669 (1995) using data from d Line 848  C Millero p.669 (1995) using data from d
848       &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +       &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +
849       &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *       &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *
850       &          dlogtk + 0.053105 _d 0*sqrts*tk)       &          dlogtk + 0.053105 _d 0*sqrts*tk)
851               if (kLevel.gt.1) then
852    C Mick and Karsten - Dec 04
853    C ADDING pressure dependence based on Millero (1995), p675
854    C with additional info from CO2sys documentation (E. Lewis and
855    C D. Wallace, 1998 - see endnotes for commentary on Millero, 95)
856                bigR = 83.145
857                dv = -29.48 + 0.1622*t + 2.608d-3*t*t
858                dk = -2.84d-3
859                pfactor = - (dv/(bigR*tk))*pressc
860         &                + (0.5*dk/(bigR*tk))*pressc*pressc
861                akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)
862               endif
863  C------------------------------------------------------------------------  C------------------------------------------------------------------------
864  C k1p = [H][H2PO4]/[H3PO4]  C k1p = [H][H2PO4]/[H3PO4]
865  C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)  C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
# Line 847  C Morris & Riley (1966) Line 920  C Morris & Riley (1966)
920  C Riley (1965)  C Riley (1965)
921             ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0             ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
922  C------------------------------------------------------------------------  C------------------------------------------------------------------------
          else  
 c add Bennington  
             fugf(i,j,bi,bj)=0. _d 0  
             ff(i,j,bi,bj)=0. _d 0  
             ak0(i,j,bi,bj)= 0. _d 0  
             ak1(i,j,bi,bj)= 0. _d 0  
             ak2(i,j,bi,bj)= 0. _d 0  
             akb(i,j,bi,bj)= 0. _d 0  
             ak1p(i,j,bi,bj) = 0. _d 0  
             ak2p(i,j,bi,bj) = 0. _d 0  
             ak3p(i,j,bi,bj) = 0. _d 0  
             aksi(i,j,bi,bj) = 0. _d 0  
             akw(i,j,bi,bj) = 0. _d 0  
             aks(i,j,bi,bj)= 0. _d 0  
             akf(i,j,bi,bj)= 0. _d 0  
             bt(i,j,bi,bj) = 0. _d 0  
             st(i,j,bi,bj) = 0. _d 0  
             ft(i,j,bi,bj) = 0. _d 0  
          endif  
          end do  
         end do  
   
         return  
         end  
   
 c=================================================================  
 c *******************************************************************  
 c=================================================================  
 CStartOfInterFace  
       SUBROUTINE CARBON_COEFFS_PRESSURE_DEP(  
      I                   ttemp,stemp,  
      I                   bi,bj,iMin,iMax,jMin,jMax,  
      I                   Klevel,myThid)  
 C  
 C     /==========================================================\  
 C     | SUBROUTINE CARBON_COEFFS                                 |  
 C     | determine coefficients for surface carbon chemistry      |  
 C     | adapted from OCMIP2:  SUBROUTINE CO2CALC                 |  
 C     | mick follows, oct 1999                                   |  
 c     | minor changes to tidy, swd aug 2002                      |  
 c     | MODIFIED FOR PRESSURE DEPENDENCE                         |  
 c     | Karsten Friis and Mick Follows 2004                      |  
 C     \==========================================================/  
 C INPUT  
 C       diclocal = total inorganic carbon (mol/m^3)  
 C             where 1 T = 1 metric ton = 1000 kg  
 C       ta  = total alkalinity (eq/m^3)  
 C       pt  = inorganic phosphate (mol/^3)  
 C       sit = inorganic silicate (mol/^3)  
 C       t   = temperature (degrees C)  
 C       s   = salinity (PSU)  
 C OUTPUT  
 C IMPORTANT: Some words about units - (JCO, 4/4/1999)  
 c     - Models carry tracers in mol/m^3 (on a per volume basis)  
 c     - Conversely, this routine, which was written by observationalists  
 c       (C. Sabine and R. Key), passes input arguments in umol/kg  
 c       (i.e., on a per mass basis)  
 c     - I have changed things slightly so that input arguments are in mol/m^3,  
 c     - Thus, all input concentrations (diclocal, ta, pt, and st) should be  
 c       given in mol/m^3; output arguments "co2star" and "dco2star"  
 c       are likewise be in mol/m^3.  
 c  
 c  
 c NOW INCLUDES:  
 c PRESSURE DEPENDENCE of K1, K2, solubility product of calcite  
 c based on Takahashi, GEOSECS Atlantic Report, Vol. 1 (1981)  
 c  
 C--------------------------------------------------------------------------  
         IMPLICIT NONE  
 C     == GLobal variables ==  
 #include "SIZE.h"  
 #include "DYNVARS.h"  
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
 #include "GRID.h"  
 #include "FFIELDS.h"  
 #include "DARWIN_FLUX.h"  
 C     == Routine arguments ==  
 C ttemp and stemp are local theta and salt arrays  
 C dont really need to pass T and S in, could use theta, salt in  
 C common block in DYNVARS.h, but this way keeps subroutine more  
 C general  
         _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
         _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
         INTEGER bi,bj,iMin,iMax,jMin,jMax  
 c K is depth index  
         INTEGER Klevel  
         INTEGER myThid  
 CEndOfInterface  
   
   
   
 C LOCAL VARIABLES  
         _RL  t  
         _RL  s  
         _RL  ta  
         _RL  pt  
         _RL  sit  
         _RL  tk  
         _RL  tk100  
         _RL  tk1002  
         _RL  dlogtk  
         _RL  sqrtis  
         _RL  sqrts  
         _RL  s15  
         _RL  scl  
         _RL  x1  
         _RL  x2  
         _RL  s2  
         _RL  xacc  
         _RL  invtk  
         _RL  is  
         _RL  is2  
         INTEGER i  
         INTEGER j  
         INTEGER k  
         _RL  bdepth  
         _RL  cdepth  
         _RL  pressc  
         _RL  Ksp_T_Calc  
         _RL  xvalue  
         _RL  zdum  
         _RL  tmpa1  
         _RL  tmpa2  
         _RL  tmpa3  
         _RL  logKspc  
         _RL  dv  
         _RL  dk  
         _RL  pfactor  
         _RL  bigR  
   
 C.....................................................................  
 C OCMIP note:  
 C Calculate all constants needed to convert between various measured  
 C carbon species. References for each equation are noted in the code.  
 C Once calculated, the constants are  
 C stored and passed in the common block "const". The original version  
 C of this code was based on the code by dickson in Version 2 of  
 C "Handbook of Methods C for the Analysis of the Various Parameters of  
 C the Carbon Dioxide System in Seawater", DOE, 1994 (SOP No. 3, p25-26).  
 C....................................................................  
   
 c determine pressure (bar) from depth  
 c 1 BAR at z=0m (atmos pressure)  
 c use UPPER surface of cell so top layer pressure = 0 bar  
 c for surface exchange coeffs  
   
 cmick..............................  
 c        write(6,*)'Klevel ',klevel  
   
         bdepth = 0.0d0  
         cdepth = 0.0d0  
         pressc = 1.0d0  
         do k = 1,Klevel  
             cdepth = bdepth + 0.5d0*drF(k)  
             bdepth = bdepth + drF(k)  
             pressc = 1.0d0 + 0.1d0*cdepth  
         end do  
 cmick...................................................  
 c        write(6,*)'depth,pressc ',cdepth,pressc  
 cmick....................................................  
   
   
   
         do i=imin,imax  
          do j=jmin,jmax  
           if (hFacC(i,j,Klevel,bi,bj).gt.0.d0) then  
            t = ttemp(i,j,Klevel,bi,bj)  
            s = max(4. _d 0, stemp(i,j,Klevel,bi,bj))  
 C terms used more than once  
            tk = 273.15 + t  
            tk100 = tk/100.0  
            tk1002=tk100*tk100  
            invtk=1.0/tk  
            dlogtk=log(tk)  
            is=19.924*s/(1000.-1.005*s)  
            is2=is*is  
            sqrtis=sqrt(is)  
            s2=s*s  
            sqrts=sqrt(s)  
            s15=s**1.5  
            scl=s/1.80655  
   
 C------------------------------------------------------------------------  
 C f = k0(1-pH2O)*correction term for non-ideality  
 C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)  
            ff(i,j,bi,bj) = exp(-162.8301 + 218.2968/tk100  +  
      &          90.9241*log(tk100) - 1.47696*tk1002 +  
      &          s * (.025695 - .025225*tk100 +  
      &          0.0049867*tk1002))  
 C------------------------------------------------------------------------  
 C K0 from Weiss 1974  
            ak0(i,j,bi,bj) = exp(93.4517/tk100 - 60.2409 +  
      &        23.3585 * log(tk100) +  
      &        s * (0.023517 - 0.023656*tk100 +  
      &        0.0047036*tk1002))  
 C------------------------------------------------------------------------  
 C k1 = [H][HCO3]/[H2CO3]  
 C k2 = [H][CO3]/[HCO3]  
 C Millero p.664 (1995) using Mehrbach et al. data on seawater scale  
            ak1(i,j,bi,bj)=10**(-1*(3670.7*invtk -  
      &          62.008 + 9.7944*dlogtk -  
      &          0.0118 * s + 0.000116*s2))  
            ak2(i,j,bi,bj)=10**(-1*(1394.7*invtk + 4.777 -  
      &          0.0184*s + 0.000118*s2))  
 C NOW PRESSURE DEPENDENCE:  
 c Following Takahashi (1981) GEOSECS report - quoting Culberson and  
 c Pytkowicz (1968)  
 c pressc = pressure in bars  
            ak1(i,j,bi,bj) = ak1(i,j,bi,bj)*  
      &             exp( (24.2-0.085*t)*(pressc-1.0)/(83.143*tk) )  
 c FIRST GO FOR K2: According to GEOSECS (1982) report  
 c          ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*  
 c    &             exp( (26.4-0.040*t)*(pressc-1.0)/(83.143*tk) )  
 c SECOND GO FOR K2: corrected coeff according to CO2sys documentation  
 c                   E. Lewis and D. Wallace (1998) ORNL/CDIAC-105  
            ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*  
      &             exp( (16.4-0.040*t)*(pressc-1.0)/(83.143*tk) )  
 C------------------------------------------------------------------------  
 C kb = [H][BO2]/[HBO2]  
 C Millero p.669 (1995) using data from dickson (1990)  
            akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s +  
      &          1.728*s15 - 0.0996*s2)*invtk +  
      &          (148.0248 + 137.1942*sqrts + 1.62142*s) +  
      &          (-24.4344 - 25.085*sqrts - 0.2474*s) *  
      &          dlogtk + 0.053105*sqrts*tk)  
 C Mick and Karsten - Dec 04  
 C ADDING pressure dependence based on Millero (1995), p675  
 C with additional info from CO2sys documentation (E. Lewis and  
 C D. Wallace, 1998 - see endnotes for commentary on Millero, 95)  
            bigR = 83.145  
            dv = -29.48 + 0.1622*t + 2.608d-3*t*t  
            dk = -2.84d-3  
            pfactor = - (dv/(bigR*tk))*pressc  
      &               + (0.5*dk/(bigR*tk))*pressc*pressc  
            akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)  
 C------------------------------------------------------------------------  
 C k1p = [H][H2PO4]/[H3PO4]  
 C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)  
            ak1p(i,j,bi,bj) = exp(-4576.752*invtk + 115.525 -  
      &          18.453*dlogtk +  
      &          (-106.736*invtk + 0.69171)*sqrts +  
      &          (-0.65643*invtk - 0.01844)*s)  
 C------------------------------------------------------------------------  
 C k2p = [H][HPO4]/[H2PO4]  
 C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))  
            ak2p(i,j,bi,bj) = exp(-8814.715*invtk + 172.0883 -  
      &          27.927*dlogtk +  
      &          (-160.340*invtk + 1.3566) * sqrts +  
      &          (0.37335*invtk - 0.05778) * s)  
 C------------------------------------------------------------------------  
 C k3p = [H][PO4]/[HPO4]  
 C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)  
            ak3p(i,j,bi,bj) = exp(-3070.75*invtk - 18.141 +  
      &          (17.27039*invtk + 2.81197) *  
      &          sqrts + (-44.99486*invtk - 0.09984) * s)  
 C------------------------------------------------------------------------  
 C ksi = [H][SiO(OH)3]/[Si(OH)4]  
 C Millero p.671 (1995) using data from Yao and Millero (1995)  
            aksi(i,j,bi,bj) = exp(-8904.2*invtk + 117.385 -  
      &          19.334*dlogtk +  
      &          (-458.79*invtk + 3.5913) * sqrtis +  
      &          (188.74*invtk - 1.5998) * is +  
      &          (-12.1652*invtk + 0.07871) * is2 +  
      &          log(1.0-0.001005*s))  
 C------------------------------------------------------------------------  
 C kw = [H][OH]  
 C Millero p.670 (1995) using composite data  
            akw(i,j,bi,bj) = exp(-13847.26*invtk + 148.9652 -  
      &          23.6521*dlogtk +  
      &          (118.67*invtk - 5.977 + 1.0495 * dlogtk) *  
      &          sqrts - 0.01615 * s)  
 C------------------------------------------------------------------------  
 C ks = [H][SO4]/[HSO4]  
 C dickson (1990, J. chem. Thermodynamics 22, 113)  
            aks(i,j,bi,bj)=exp(-4276.1*invtk + 141.328 -  
      &          23.093*dlogtk +  
      &          (-13856*invtk + 324.57 - 47.986*dlogtk)*sqrtis +  
      &          (35474*invtk - 771.54 + 114.723*dlogtk)*is -  
      &          2698*invtk*is**1.5 + 1776*invtk*is2 +  
      &          log(1.0 - 0.001005*s))  
 C------------------------------------------------------------------------  
 C kf = [H][F]/[HF]  
 C dickson and Riley (1979) -- change pH scale to total  
            akf(i,j,bi,bj)=exp(1590.2*invtk - 12.641 + 1.525*sqrtis +  
      &          log(1.0 - 0.001005*s) +  
      &          log(1.0 + (0.1400/96.062)*(scl)/aks(i,j,bi,bj)))  
 C------------------------------------------------------------------------  
 C Calculate concentrations for borate, sulfate, and fluoride  
 C Uppstrom (1974)  
            bt(i,j,bi,bj) = 0.000232 * scl/10.811  
 C Morris & Riley (1966)  
            st(i,j,bi,bj) = 0.14 * scl/96.062  
 C Riley (1965)  
            ft(i,j,bi,bj) = 0.000067 * scl/18.9984  
 C------------------------------------------------------------------------  
923  C solubility product for calcite  C solubility product for calcite
924  C  C
925  c Following Takahashi (1982) GEOSECS handbook  c Following Takahashi (1982) GEOSECS handbook
# Line 1186  c alternative pressure dependence from I Line 963  c alternative pressure dependence from I
963    
964             Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue)             Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue)
965    
   
   
   
966  C------------------------------------------------------------------------  C------------------------------------------------------------------------
967           else           else
968              ff(i,j,bi,bj)=0.d0  c add Bennington
969              ak0(i,j,bi,bj)= 0.d0              fugf(i,j,bi,bj)=0. _d 0
970              ak1(i,j,bi,bj)= 0.d0              ff(i,j,bi,bj)=0. _d 0
971              ak2(i,j,bi,bj)= 0.d0              ak0(i,j,bi,bj)= 0. _d 0
972              akb(i,j,bi,bj)= 0.d0              ak1(i,j,bi,bj)= 0. _d 0
973              ak1p(i,j,bi,bj) = 0.d0              ak2(i,j,bi,bj)= 0. _d 0
974              ak2p(i,j,bi,bj) = 0.d0              akb(i,j,bi,bj)= 0. _d 0
975              ak3p(i,j,bi,bj) = 0.d0              ak1p(i,j,bi,bj) = 0. _d 0
976              aksi(i,j,bi,bj) = 0.d0              ak2p(i,j,bi,bj) = 0. _d 0
977              akw(i,j,bi,bj) = 0.d0              ak3p(i,j,bi,bj) = 0. _d 0
978              aks(i,j,bi,bj)= 0.d0              aksi(i,j,bi,bj) = 0. _d 0
979              akf(i,j,bi,bj)= 0.d0              akw(i,j,bi,bj) = 0. _d 0
980              bt(i,j,bi,bj) = 0.d0              aks(i,j,bi,bj)= 0. _d 0
981              st(i,j,bi,bj) = 0.d0              akf(i,j,bi,bj)= 0. _d 0
982              ft(i,j,bi,bj) = 0.d0              bt(i,j,bi,bj) = 0. _d 0
983                st(i,j,bi,bj) = 0. _d 0
984                ft(i,j,bi,bj) = 0. _d 0
985              Ksp_TP_Calc(i,j,bi,bj) = 0.d0              Ksp_TP_Calc(i,j,bi,bj) = 0.d0
986           endif           endif
987           end do           end do
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22