--- MITgcm/pkg/seaice/seaice_jfnk.F	2013/05/30 14:07:19	1.23
+++ MITgcm/pkg/seaice/seaice_jfnk.F	2014/12/01 12:31:36	1.28
@@ -1,7 +1,10 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.23 2013/05/30 14:07:19 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.28 2014/12/01 12:31:36 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
+#ifdef ALLOW_AUTODIFF
+# include "AUTODIFF_OPTIONS.h"
+#endif
 
 C--  File seaice_jfnk.F: seaice jfnk dynamical solver S/R:
 C--   Contents
@@ -76,8 +79,8 @@
       _RL     JFNKgamma_lin
       _RL     FGMRESeps
       _RL     JFNKtol
-C     Adams-Bashforth extrapolation factors
-      _RL abFac, abAlpha
+C     backward differences extrapolation factors
+      _RL bdfFac, bdfAlpha
 C
       _RL     recip_deltaT
       LOGICAL JFNKconverged, krylovConverged
@@ -87,7 +90,7 @@
 C     u/vIceRes :: residual of sea-ice momentum equations
       _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
-C     extra time level required for Adams-Bashforth-2 time stepping
+C     extra time level required for backward difference time stepping
       _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 C     du/vIce   :: ice velocity increment to be added to u/vIce
@@ -96,6 +99,7 @@
 C     precomputed (= constant per Newton iteration) versions of
 C     zeta, eta, and DWATN, press
       _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
@@ -118,16 +122,16 @@
      &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
      &     iOutFGMRES=1
 
-C     Adams-Bashforth extrapolation factors
-      abFac = 0. _d 0
-      IF ( SEAICEuseAB2 ) THEN
-       IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartAB.EQ.0 ) THEN
-        abFac = 0. _d 0
+C     backward difference extrapolation factors
+      bdfFac = 0. _d 0
+      IF ( SEAICEuseBDF2 ) THEN
+       IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
+        bdfFac = 0. _d 0
        ELSE
-        abFac = 0.5 _d 0 + SEAICE_abEps
+        bdfFac = 0.5 _d 0
        ENDIF
       ENDIF
-      abAlpha = 1. _d 0 + abFac 
+      bdfAlpha = 1. _d 0 + bdfFac 
 
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
@@ -142,18 +146,19 @@
 C     cycle ice velocities
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
-          duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * abAlpha 
-     &         + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * abFac
-          dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * abAlpha 
-     &         + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * abFac
+          duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha 
+     &         + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac
+          dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha 
+     &         + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac
           uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
           vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
          ENDDO
         ENDDO
-        IF ( .NOT.SEAICEuseIMEX ) THEN
+C     As long as IMEX is not properly implemented leave this commented out
+CML        IF ( .NOT.SEAICEuseIMEX ) THEN
 C     Compute things that do no change during the Newton iteration:
 C     sea-surface tilt and wind stress:
-C     FORCEX/Y0 - mass*(abA*u/vIceNm1+abB*(u/vIceNm1-u/vIceNm2))/deltaT
+C     FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj)
@@ -162,7 +167,7 @@
      &         + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
          ENDDO
         ENDDO
-        ENDIF
+CML        ENDIF
        ENDDO
       ENDDO
 C     Start nonlinear Newton iteration: outer loop iteration
@@ -183,6 +188,7 @@
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)
+           zetaZPre(I,J,bi,bj)= zetaZ(I,J,bi,bj)
             etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)
            etaZPre(I,J,bi,bj) =  etaZ(I,J,bi,bj)
            dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
@@ -239,7 +245,7 @@
           IF ( SOLV_MAX_ITERS .GT. 0 )
      &         CALL SEAICE_PRECONDITIONER(
      U         duIce, dvIce,
-     I         zetaPre, etaPre, etaZpre, dwatPre,
+     I         zetaPre, etaPre, etaZpre, zetaZpre, dwatPre,
      I         newtonIter, krylovIter, myTime, myIter, myThid )
          ELSEIF (iCode.GE.2) THEN
 C     Compute Jacobian times vector