--- MITgcm/pkg/seaice/seaice_jfnk.F	2014/10/20 03:20:57	1.27
+++ MITgcm/pkg/seaice/seaice_jfnk.F	2016/01/27 14:03:34	1.29
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.27 2014/10/20 03:20:57 gforget Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.29 2016/01/27 14:03:34 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
@@ -63,12 +63,18 @@
 
 C     !LOCAL VARIABLES:
 C     === Local variables ===
-C     i,j,bi,bj :: loop indices
-      INTEGER i,j,bi,bj
+C     i,j,k,bi,bj :: loop indices
+      INTEGER i,j,k,bi,bj
 C     loop indices
       INTEGER newtonIter
       INTEGER krylovIter, krylovFails
       INTEGER totalKrylovItersLoc, totalNewtonItersLoc
+C     FGMRES parameters
+C     im      :: size of Krylov space
+C     ifgmres :: interation counter
+      INTEGER im
+      PARAMETER ( im = 50 )
+      INTEGER ifgmres
 C     FGMRES flag that determines amount of output messages of fgmres
       INTEGER iOutFGMRES
 C     FGMRES flag that indicates what fgmres wants us to do next
@@ -99,9 +105,14 @@
 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)
+C     work arrays
+      _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy)
+      _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy)
+      _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy)
 CEOP
 
 C     Initialise
@@ -187,6 +198,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)
@@ -225,16 +237,36 @@
 
         krylovConverged = .FALSE.
         FGMRESeps = JFNKgamma_lin * JFNKresidual
+C     map first guess sol; it is zero because the solution is a correction
+       CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid)
+C     map rhs and change its sign because we are solving J*u = -F
+        CALL SEAICE_MAP2VEC(nVec,-uIceRes,-vIceRes,rhs,.TRUE.,myThid)
         DO WHILE ( .NOT.krylovConverged )
 C     solution vector sol = du/vIce
 C     residual vector (rhs) Fu = u/vIceRes
 C     output work vectors wk1, -> input work vector wk2
 
-         CALL SEAICE_FGMRES_DRIVER(
-     I        uIceRes, vIceRes,
-     U        duIce, dvIce, iCode,
-     I        FGMRESeps, iOutFGMRES,
-     I        newtonIter, krylovIter, myTime, myIter, myThid )
+C     map preconditioner results or Jacobian times vector,
+C     stored in du/vIce to wk2, for iCode=0, wk2 is set to zero,
+C     because du/vIce = 0
+         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid)
+C
+         CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,
+     U        vv,w,wk1,wk2,
+     I        FGMRESeps,SEAICEkrylovIterMax,iOutFGMRES,
+     U        iCode,
+     I        myThid)
+C     
+         IF ( iCode .EQ. 0 ) THEN
+C     map sol(ution) vector to du/vIce
+          CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid)
+         ELSE
+C     map work vector to du/vIce to either compute a preconditioner
+C     solution (wk1=rhs) or a Jacobian times wk1
+          CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid)
+         ENDIF
+C     Fill overlaps in updated fields
+         CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)
 C     FGMRES returns iCode either asking for an new preconditioned vector
 C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
 C     iteration
@@ -243,7 +275,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