--- MITgcm/pkg/dic/dic_surfforcing.F	2003/06/25 21:00:36	1.1
+++ MITgcm/pkg/dic/dic_surfforcing.F	2003/10/09 04:19:19	1.4
@@ -1,4 +1,4 @@
-#include "CPP_OPTIONS.h"
+#include "DIC_OPTIONS.h"
 #include "PTRACERS_OPTIONS.h"
 #include "GCHEM_OPTIONS.h"
 
@@ -34,13 +34,9 @@
       INTEGER iMin,iMax,jMin,jMax, bi, bj
 
 #ifdef ALLOW_PTRACERS
-#ifdef DIC_ABIOTIC
 C     == Local variables ==
-       INTEGER I,J, kLev
+       INTEGER I,J, kLev, it
 C Number of iterations for pCO2 solvers...
-      INTEGER inewtonmax
-      INTEGER ibrackmax
-      INTEGER donewt
 C Solubility relation coefficients
       _RL SchmidtNoDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL pCO2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
@@ -50,7 +46,6 @@
       _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
-      _RL FluxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
@@ -81,7 +76,7 @@
              surfphos(i,j)  = 5.1225e-4 * maskC(i,j,kLev,bi,bj)
 #endif
 C FOR NON-INTERACTIVE Si
-             surfsi(i,j)   = 7.6838e-3 * maskC(i,j,kLev,bi,bj)
+             surfsi(i,j)   = SILICA(i,j,bi,bj) * maskC(i,j,kLev,bi,bj)
           ENDDO
          ENDDO
 
@@ -90,59 +85,13 @@
      I                       bi,bj,iMin,iMax,jMin,jMax)
 C====================================================================
 
-#define PH_APPROX
-c set number of iterations for [H+] solvers
-#ifdef PH_APPROX
-       inewtonmax = 1
-#else
-       inewtonmax = 10
-#endif
-       ibrackmax = 30
-C determine pCO2 in surface ocean
-C set guess of pH for first step here
-C IF first step THEN use bracket-bisection for first step, 
-C and determine carbon coefficients for safety
-C ELSE use newton-raphson with previous H+(x,y) as first guess
-
-       donewt=1
-
-c for first few timesteps
-       IF(myIter .le. (nIter0+inewtonmax) )then
-          donewt=0
-          DO j=1-OLy,sNy+OLy
-           DO i=1-OLx,sNx+OLx
-                  pH(i,j,bi,bj) = 8.0
-           ENDDO
-          ENDDO
-#ifdef PH_APPROX
-          print*,'QQ: pCO2 approximation method'
-c first approxmation
-       DO j=1-OLy,sNy+OLy
-        DO i=1-OLx,sNx+OLx
-          CALL CALC_PCO2_APPROX(
-     I        theta(i,j,kLev,bi,bj),salt(i,j,kLev,bi,bj),
-     I        PTR_CO2(i,j,kLev), surfphos(i,j),
-     I        surfsi(i,j),surfalk(i,j),
-     I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),
-     I        ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
-     I        aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
-     I        aksi(i,j,bi,bj),akf(i,j,bi,bj),ff(i,j,bi,bj),
-     I        bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
-     U        pH(i,j,bi,bj),pCO2(i,j,bi,bj) )
-        ENDDO
-       ENDDO
-#else
-          print*,'QQ: pCO2 full method'
-#endif
-       ENDIF
-
-
 c pCO2 solver...
+C$TAF LOOP = parallel
        DO j=1-OLy,sNy+OLy
+C$TAF LOOP = parallel
         DO i=1-OLx,sNx+OLx
 
           IF(maskC(i,j,kLev,bi,bj) .NE. 0.)THEN
-#ifdef PH_APPROX
             CALL CALC_PCO2_APPROX(
      I        theta(i,j,kLev,bi,bj),salt(i,j,kLev,bi,bj),
      I        PTR_CO2(i,j,kLev), surfphos(i,j),
@@ -153,18 +102,6 @@
      I        aksi(i,j,bi,bj),akf(i,j,bi,bj),ff(i,j,bi,bj),
      I        bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
      U        pH(i,j,bi,bj),pCO2(i,j,bi,bj) )
-#else
-            CALL CALC_PCO2(donewt,inewtonmax,ibrackmax,
-     I        theta(i,j,kLev,bi,bj),salt(i,j,kLev,bi,bj),
-     I        PTR_CO2(i,j,kLev), surfphos(i,j),
-     I        surfsi(i,j),surfalk(i,j),
-     I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),
-     I        ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
-     I        aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
-     I        aksi(i,j,bi,bj),akf(i,j,bi,bj),ff(i,j,bi,bj),
-     I        bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
-     U        pH(i,j,bi,bj),pCO2(i,j,bi,bj) )
-#endif
           ELSE
              pCO2(i,j,bi,bj)=0. _d 0
           END IF
@@ -202,16 +139,16 @@
 C Flux = Vp * ([CO2sat] - [CO2])
 C CO2sat = K0*pCO2atmos*P/P0
 C Converting pCO2 to [CO2] using ff, as in CALC_PCO2
-              FluxCO2(i,j) = 
+              FluxCO2(i,j,bi,bj) = 
      &         maskC(i,j,kLev,bi,bj)*Kwexch(i,j)*( 
      &         ak0(i,j,bi,bj)*pCO2sat(i,j) - 
      &         ff(i,j,bi,bj)*pCO2(i,j,bi,bj) 
      &         ) 
             ELSE
-               FluxCO2(i,j) = 0.
+               FluxCO2(i,j,bi,bj) = 0.
             ENDIF
 C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
-            FluxCO2(i,j) = FluxCO2(i,j)/permil
+            FluxCO2(i,j,bi,bj) = FluxCO2(i,j,bi,bj)/permil
 
             IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
 c calculate virtual flux
@@ -236,13 +173,12 @@
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            GDC(i,j)= maskC(i,j,kLev,bi,bj)*(
-     &                    FluxCO2(i,j)*recip_drF(kLev)
+     &                    FluxCO2(i,j,bi,bj)*recip_drF(kLev)
      &                    + VirtualFlux(i,j)
      &                                              )
           ENDDO
          ENDDO
 
 #endif
-#endif
         RETURN
         END