--- MITgcm/pkg/dic/fe_chem.F	2003/10/09 04:19:19	1.3
+++ MITgcm/pkg/dic/fe_chem.F	2006/11/28 21:16:03	1.9
@@ -1,7 +1,9 @@
+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 $Name:  $
+
 #include "DIC_OPTIONS.h"
 #include "GCHEM_OPTIONS.h"
 
-#ifdef ALLOW_FE
 CStartOfInterFace
       SUBROUTINE Fe_CHEM(
      I           bi,bj,iMin,iMax,jMin,jMax,
@@ -21,7 +23,6 @@
 #include "PARAMS.h"
 #include "GRID.h"
 #include "DIC_BIOTIC.h"
-#include "PTRACERS.h"
 
 C     == Routine arguments ==                
 C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
@@ -31,8 +32,12 @@
       _RL  fe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       INTEGER bi,bj,iMin,iMax,jMin,jMax
       INTEGER myIter,myThid
+#ifdef AD_SAFE
+      _RL thx, thy, theps
+#endif
 CEndOfInterface
 
+#ifdef ALLOW_FE
 
       INTEGER I,J,K
       _RL  lig, FeL
@@ -51,8 +56,8 @@
 C ligand balance in surface layer
 C in surface layer
 
-       DO j=jMin,jMax
-        DO i=iMin,iMax
+       DO j=jmin,jmax
+        DO i=imin,imax
          DO k=1,nR
           IF(hFacC(i,j,k,bi,bj) .gt. 0.0)THEN
 
@@ -66,11 +71,25 @@
  
               FeL = ligand_tot-lig
               freefe(i,j,k,bi,bj) = fe (i,j,k,bi,bj)-FeL
+#ifdef MINFE
+#ifdef AD_SAFE
+              thx=freefe(i,j,k,bi,bj)
+              thy=freefemax
+              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
+
+#else
+              freefe(i,j,k,bi,bj) = min(freefe(i,j,k,bi,bj),freefemax)
+#endif
+              fe(i,j,k,bi,bj) = FeL+freefe(i,j,k,bi,bj)
+#endif
           END IF  
          ENDDO
         ENDDO
        ENDDO    
 c
+#endif 
         RETURN
         END                      
-#endif