--- MITgcm/pkg/dic/fe_chem.F	2006/11/28 21:16:03	1.9
+++ MITgcm/pkg/dic/fe_chem.F	2008/04/07 20:31:16	1.13
@@ -1,8 +1,7 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/dic/fe_chem.F,v 1.9 2006/11/28 21:16:03 stephd Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/dic/fe_chem.F,v 1.13 2008/04/07 20:31:16 dfer Exp $
 C $Name:  $
 
 #include "DIC_OPTIONS.h"
-#include "GCHEM_OPTIONS.h"
 
 CStartOfInterFace
       SUBROUTINE Fe_CHEM(
@@ -22,7 +21,7 @@
 #include "EEPARAMS.h"
 #include "PARAMS.h"
 #include "GRID.h"
-#include "DIC_BIOTIC.h"
+#include "DIC_VARS.h"
 
 C     == Routine arguments ==                
 C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
@@ -41,6 +40,7 @@
 
       INTEGER I,J,K
       _RL  lig, FeL
+      _RL  tmpfe
 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
@@ -58,27 +58,37 @@
 
        DO j=jmin,jmax
         DO i=imin,imax
-         DO k=1,nR
-          IF(hFacC(i,j,k,bi,bj) .gt. 0.0)THEN
+         DO k=1,Nr
+          IF (hFacC(i,j,k,bi,bj) .GT. 0. _d 0) THEN
 
-C   Ligand,FeL,Fe calculation
+#ifdef DIC_NO_NEG
+              tmpfe=max(0. _d0 , fe (i,j,k,bi,bj))
+#else
+              tmpfe=fe (i,j,k,bi,bj)
+#endif
             
-              lig=(-ligand_stab*fe (i,j,k,bi,bj)+
-     &              ligand_stab*ligand_tot-1
-     &             +((ligand_stab*fe (i,j,k,bi,bj)
-     &                -ligand_stab*ligand_tot+1)**2+4
-     &                *ligand_stab*ligand_tot)**0.5)/(2*ligand_stab)
+C   Ligand,FeL,Fe calculation
+              lig=(-ligand_stab*tmpfe +
+     &              ligand_stab*ligand_tot-1. _d 0
+     &             +((ligand_stab*tmpfe 
+     &                -ligand_stab*ligand_tot+1. _d 0)**2
+     &               +4. _d 0*ligand_stab*ligand_tot)**0.5 _d 0
+     &            )/(2. _d 0*ligand_stab)
  
               FeL = ligand_tot-lig
-              freefe(i,j,k,bi,bj) = fe (i,j,k,bi,bj)-FeL
+              if (tmpfe.eq.0. _d 0) then
+                freefe(i,j,k,bi,bj) = tmpfe -FeL
+              else
+                freefe(i,j,k,bi,bj) = 0. _d 0
+              endif
 #ifdef MINFE
 #ifdef AD_SAFE
               thx=freefe(i,j,k,bi,bj)
               thy=freefemax
-              theps=1.d-8
+              theps=1. _d -8
               freefe(i,j,k,bi,bj) =
-     &                 ( 1.d0 - tanh((thx-thy)/theps) ) * thx/2 +
-     &                 ( 1.d0 + tanh((thx-thy)/theps) ) * thy/2
+     &                 ( 1. _d 0 - tanh((thx-thy)/theps) ) * thx/2.+
+     &                 ( 1. _d 0 + tanh((thx-thy)/theps) ) * thy/2.
 
 #else
               freefe(i,j,k,bi,bj) = min(freefe(i,j,k,bi,bj),freefemax)
@@ -92,4 +102,4 @@
 c
 #endif 
         RETURN
-        END                      
+        END