--- MITgcm/pkg/seaice/seaice_jfnk.F	2013/04/04 07:02:51	1.21
+++ MITgcm/pkg/seaice/seaice_jfnk.F	2014/02/07 14:27:21	1.25
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.21 2013/04/04 07:02:51 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.25 2014/02/07 14:27:21 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
@@ -73,10 +73,12 @@
       _RL     JFNKresidual
       _RL     JFNKresidualKm1
 C     parameters to compute convergence criterion
-      _RL     phi_e, alp_e, JFNKgamma_lin
+      _RL     JFNKgamma_lin
       _RL     FGMRESeps
       _RL     JFNKtol
-
+C     backward differences extrapolation factors
+      _RL bdfFac, bdfAlpha
+C
       _RL     recip_deltaT
       LOGICAL JFNKconverged, krylovConverged
       LOGICAL writeNow
@@ -85,6 +87,9 @@
 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 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
       _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
@@ -113,6 +118,17 @@
      &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
      &     iOutFGMRES=1
 
+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
+        bdfFac = 0.5 _d 0
+       ENDIF
+      ENDIF
+      bdfAlpha = 1. _d 0 + bdfFac 
+
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO J=1-OLy,sNy+OLy
@@ -121,21 +137,32 @@
           vIceRes(I,J,bi,bj) = 0. _d 0
           duIce  (I,J,bi,bj) = 0. _d 0
           dvIce  (I,J,bi,bj) = 0. _d 0
+         ENDDO
+        ENDDO
+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) * 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     Compute things that do no change during the Newton iteration:
 C     sea-surface tilt and wind stress:
-C     FORCEX/Y0 - mass*(u/vIceNm1)/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)
-     &         + seaiceMassU(I,J,bi,bj)*uIceNm1(I,J,bi,bj)*recip_deltaT
+     &         + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT
           FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
-     &         + seaiceMassV(I,J,bi,bj)*vIceNm1(I,J,bi,bj)*recip_deltaT
+     &         + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
          ENDDO
         ENDDO
+        ENDIF
        ENDDO
       ENDDO
 C     Start nonlinear Newton iteration: outer loop iteration
@@ -167,10 +194,9 @@
        JFNKgamma_lin = JFNKgamma_lin_max
        IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter
      &      .AND.JFNKresidual.LT.JFNKres_t ) THEN
-C     Eisenstat, 1996, equ.(2.6)
-        phi_e = 1. _d 0
-        alp_e = 1. _d 0
-        JFNKgamma_lin = phi_e*( JFNKresidual/JFNKresidualKm1 )**alp_e
+C     Eisenstat and Walker (1996), eq.(2.6)
+        JFNKgamma_lin = SEAICE_JFNKphi
+     &       *( JFNKresidual/JFNKresidualKm1 )**SEAICE_JFNKalpha
         JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin)
         JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin)
        ENDIF