C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.17 2013/01/17 10:42:43 mlosch Exp $
C $Name: checkpoint64c $

#include "SEAICE_OPTIONS.h"

C--  File seaice_jfnk.F: seaice jfnk dynamical solver S/R:
C--   Contents
C--   o SEAICE_JFNK
C--   o SEAICE_JFNK_UPDATE

CBOP
C     !ROUTINE: SEAICE_JFNK
C     !INTERFACE:
      SUBROUTINE SEAICE_JFNK( myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SEAICE_JFNK
C     | o Ice dynamics using a Jacobian-free Newton-Krylov solver
C     |   following J.-F. Lemieux et al. Improving the numerical
C     |   convergence of viscous-plastic sea ice models with the
C     |   Jacobian-free Newton-Krylov method. J. Comp. Phys. 229,
C     |   2840-2852 (2010).
C     | o The logic follows JFs code.
C     *==========================================================*
C     | written by Martin Losch, Oct 2012
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"

#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: my Thread Id. number
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#if ( (defined SEAICE_CGRID) && \
      (defined SEAICE_ALLOW_JFNK) && \
      (defined SEAICE_ALLOW_DYNAMICS) )
C     !FUNCTIONS:
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !LOCAL VARIABLES:
C     === Local variables ===
C     i,j,bi,bj :: loop indices
      INTEGER i,j,bi,bj
C     loop indices
      INTEGER newtonIter
      INTEGER krylovIter, krylovFails
      INTEGER totalKrylovItersLoc, totalNewtonItersLoc
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
      INTEGER iCode
      _RL     JFNKresidual
      _RL     JFNKresidualKm1
C     parameters to compute convergence criterion
      _RL     phi_e, alp_e, JFNKgamma_lin
      _RL     FGMRESeps
      _RL     JFNKtol
C     
      _RL     recip_deltaT
      LOGICAL JFNKconverged, krylovConverged
      LOGICAL writeNow
      CHARACTER*(MAX_LEN_MBUF) msgBuf
C
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     vector version of the residuals
      _RL resTmp (nVec,1,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)
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 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)
CEOP

C     Initialise
      newtonIter          = 0
      krylovFails         = 0
      totalKrylovItersLoc = 0
      JFNKconverged       = .FALSE.
      JFNKtol             = 0. _d 0
      JFNKresidual        = 0. _d 0
      JFNKresidualKm1     = 0. _d 0
      FGMRESeps           = 0. _d 0
      recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn

      iOutFGMRES=0
C     with iOutFgmres=1, seaice_fgmres prints the residual at each iteration
      IF ( debugLevel.GE.debLevC .AND.
     &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
     &     iOutFGMRES=1

C     
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          uIceRes(I,J,bi,bj) = 0. _d 0
          vIceRes(I,J,bi,bj) = 0. _d 0
          duIce  (I,J,bi,bj) = 0. _d 0
          dvIce  (I,J,bi,bj) = 0. _d 0
          uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
          vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
         ENDDO
        ENDDO
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
        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    
          FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
     &         + seaiceMassV(I,J,bi,bj)*vIceNm1(I,J,bi,bj)*recip_deltaT    
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C     Start nonlinear Newton iteration: outer loop iteration
      DO WHILE ( newtonIter.LT.SEAICEnewtonIterMax .AND.
     &     .NOT.JFNKconverged )
       newtonIter = newtonIter + 1
C     Compute initial residual F(u), (includes computation of global
C     variables DWATN, zeta, and eta)
       IF ( newtonIter .EQ. 1 ) CALL SEAICE_JFNK_UPDATE( 
     I      duIce, dvIce, 
     U      uIce, vIce, JFNKresidual,
     O      uIceRes, vIceRes,
     I      newtonIter, myTime, myIter, myThid )
C     local copies of precomputed coefficients that are to stay
C     constant for the preconditioner
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           zetaPre(I,J,bi,bj) =  zeta(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)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
C     compute convergence criterion for linear preconditioned FGMRES
       JFNKgamma_lin = JFNKgamma_lin_max
       IF ( newtonIter.GT.1.AND.newtonIter.LE.100
     &      .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
        JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin)
        JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin)
       ENDIF
C     save the residual for the next iteration
       JFNKresidualKm1 = JFNKresidual
C
C     The Krylov iteration using FGMRES, the preconditioner is LSOR
C     for now. The code is adapted from SEAICE_LSR, but heavily stripped
C     down.
C     krylovIter is mapped into "its" in seaice_fgmres and is incremented
C     in that routine
       krylovIter    = 0
       iCode         = 0
C
       JFNKconverged = JFNKresidual.LT.JFNKtol
C
C     do Krylov loop only if convergence is not reached
C
       IF ( .NOT.JFNKconverged ) THEN
C
C     start Krylov iteration (FGMRES)
C
        krylovConverged = .FALSE.
        FGMRESeps = JFNKgamma_lin * JFNKresidual
        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 
C     
         CALL SEAICE_FGMRES_DRIVER(
     I        uIceRes, vIceRes, 
     U        duIce, dvIce, iCode,
     I        FGMRESeps, iOutFGMRES,
     I        newtonIter, krylovIter, myTime, myIter, 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
         IF (iCode.EQ.1) THEN
C     Call preconditioner 
          IF ( SOLV_MAX_ITERS .GT. 0 )
     &         CALL SEAICE_PRECONDITIONER( 
     U         duIce, dvIce, 
     I         zetaPre, etaPre, etaZpre, dwatPre, 
     I         newtonIter, krylovIter, myTime, myIter, myThid )
         ELSEIF (iCode.GE.2) THEN
C     Compute Jacobian times vector
          CALL SEAICE_JACVEC(
     I         uIce, vIce, uIceRes, vIceRes,
     U         duIce, dvIce,  
     I         newtonIter, krylovIter, myTime, myIter, myThid )
         ENDIF
         krylovConverged = iCode.EQ.0
C     End of Krylov iterate
        ENDDO
        totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
C     some output diagnostics
        IF ( debugLevel.GE.debLevA ) THEN
         _BEGIN_MASTER( myThid )
         totalNewtonItersLoc = 
     &        SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter
         WRITE(msgBuf,'(2A,2(1XI6),2E12.5)') 
     &        ' S/R SEAICE_JFNK: Newton iterate / total, ',
     &        'JFNKgamma_lin, initial norm = ',
     &        newtonIter, totalNewtonItersLoc,
     &        JFNKgamma_lin,JFNKresidual
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &        SQUEEZE_RIGHT, myThid )
         WRITE(msgBuf,'(3(A,I6))')
     &        ' S/R SEAICE_JFNK: Newton iterate / total = ',newtonIter, 
     &        ' / ', totalNewtonItersLoc,
     &        ', Nb. of FGMRES iterations = ', krylovIter
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &        SQUEEZE_RIGHT, myThid )
         _END_MASTER( myThid )
        ENDIF
        IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
         krylovFails = krylovFails + 1
        ENDIF
C     Set the stopping criterion for the Newton iteration and the
C     criterion for the transition from accurate to approximate FGMRES
        IF ( newtonIter .EQ. 1 ) THEN
         JFNKtol=JFNKgamma_nonlin*JFNKresidual
         IF ( JFNKres_tFac .NE. UNSET_RL )
     &        JFNKres_t = JFNKresidual * JFNKres_tFac
        ENDIF
C     Update linear solution vector and return to Newton iteration
C     Do a linesearch if necessary, and compute a new residual.
C     Note that it should be possible to do the following operations
C     at the beginning of the Newton iteration, thereby saving us from
C     the extra call of seaice_jfnk_update, but unfortunately that
C     changes the results, so we leave the stuff here for now.
        CALL SEAICE_JFNK_UPDATE( 
     I       duIce, dvIce, 
     U       uIce, vIce, JFNKresidual,
     O       uIceRes, vIceRes,
     I       newtonIter, myTime, myIter, myThid )
C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO J=1-Oly,sNy+Oly
           DO I=1-Olx,sNx+Olx
            duIce(I,J,bi,bj)= 0. _d 0
            dvIce(I,J,bi,bj)= 0. _d 0
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF
C     end of Newton iterate
      ENDDO
C
C--   Output diagnostics
C
      IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
C     Count iterations
       totalJFNKtimeSteps = totalJFNKtimeSteps + 1
       totalNewtonIters   = totalNewtonIters + newtonIter
       totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
C     Record failure
       totalKrylovFails   = totalKrylovFails + krylovFails
       IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
        totalNewtonFails = totalNewtonFails + 1 
       ENDIF
      ENDIF
C     Decide whether it is time to dump and reset the counter
      writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
     &     myTime+deltaTClock, deltaTClock) 
#ifdef ALLOW_CAL
      IF ( useCAL ) THEN
       CALL CAL_TIME2DUMP( 
     I      zeroRL, SEAICE_monFreq,  deltaTClock,
     U      writeNow,
     I      myTime+deltaTclock, myIter+1, myThid )
      ENDIF
#endif
      IF ( writeNow ) THEN
       _BEGIN_MASTER( myThid )
       WRITE(msgBuf,'(A)') 
     &' // ======================================================='
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A)') 
     &' // ======================================================='
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: time step              = ', myIter+1
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A,I10)') 
     &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A)') 
     &' // ======================================================='
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A)') ' // End JFNK statistics'
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       WRITE(msgBuf,'(A)') 
     &' // ======================================================='
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       _END_MASTER( myThid )
C     reset and start again
       totalJFNKtimeSteps = 0
       totalNewtonIters   = 0
       totalKrylovIters   = 0
       totalKrylovFails   = 0
       totalNewtonFails   = 0
      ENDIF

C     Print more debugging information
      IF ( debugLevel.GE.debLevA ) THEN
       IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A,I10)') 
     &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
     &       myIter+1
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &       SQUEEZE_RIGHT, myThid )
        _END_MASTER( myThid )
       ENDIF
       IF ( krylovFails .GT. 0 ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A,I4,A,I10)') 
     &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
     &       krylovFails, ' times in timestep ', myIter+1
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &       SQUEEZE_RIGHT, myThid )
        _END_MASTER( myThid )
       ENDIF
       _BEGIN_MASTER( myThid )
       WRITE(msgBuf,'(A,I6,A,I10)') 
     &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
     &      totalKrylovItersLoc, ' in timestep ', myIter+1
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
       _END_MASTER( myThid )
      ENDIF

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: SEAICE_JFNK_UPDATE
C     !INTERFACE:

      SUBROUTINE SEAICE_JFNK_UPDATE( 
     I     duIce, dvIce, 
     U     uIce, vIce, JFNKresidual,
     O     uIceRes, vIceRes,
     I     newtonIter, myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SEAICE_JFNK_UPDATE
C     | o Update velocities with incremental solutions of FGMRES
C     | o compute residual of updated solutions and do
C     | o linesearch:
C     |   reduce update until residual is smaller than previous
C     |   one (input)
C     *==========================================================*
C     | written by Martin Losch, Jan 2013
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: my Thread Id. number
C     newtonIter :: current iterate of Newton iteration
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER newtonIter
C     JFNKresidual :: Residual at the beginning of the FGMRES iteration,
C                     changes with newtonIter (updated)
      _RL     JFNKresidual
C     du/vIce   :: ice velocity increment to be added to u/vIce (input)
      _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     u/vIce    :: ice velocity increment to be added to u/vIce (updated)
      _RL uIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     u/vIceRes :: residual of sea-ice momentum equations (output)
      _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     !LOCAL VARIABLES:
C     === Local variables ===
C     i,j,bi,bj :: loop indices
      INTEGER i,j,bi,bj
      INTEGER l
      _RL     resLoc, facLS
      LOGICAL doLineSearch
C     nVec    :: size of the input vector(s)
C     vector version of the residuals
      INTEGER nVec
      PARAMETER ( nVec  = 2*sNx*sNy )
      _RL resTmp (nVec,1,nSx,nSy)
C
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

C     Initialise some local variables
      l = 0
      resLoc = JFNKresidual
      facLS = 1. _d 0
      doLineSearch = .TRUE.
      DO WHILE ( doLineSearch )
C     Determine, if we need more iterations
       doLineSearch = resLoc .GE. JFNKresidual 
C     Limit the maximum number of iterations arbitrarily to four
       doLineSearch = doLineSearch .AND. l .LE. 4 
C     For the first iteration du/vIce = 0 and there will be no
C     improvement of the residual possible, so we do only the first
C     iteration
       IF ( newtonIter .EQ. 1 ) doLineSearch = .FALSE.
C     Only start a linesearch after some Newton iterations
       IF ( newtonIter .LE. SEAICE_JFNK_lsIter ) doLineSearch = .FALSE.
C     Increment counter
       l = l + 1
C     Create update
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
           uIce(I,J,bi,bj) = uIce(I,J,bi,bj)+facLS*duIce(I,J,bi,bj)
           vIce(I,J,bi,bj) = vIce(I,J,bi,bj)+facLS*dvIce(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
C     Compute current residual F(u), (includes re-computation of global
C     variables DWATN, zeta, and eta, i.e. they are different after this)
       CALL SEAICE_CALC_RESIDUAL( 
     I      uIce, vIce, 
     O      uIceRes, vIceRes, 
     I      newtonIter, 0, myTime, myIter, myThid )
C     Important: Compute the norm of the residual using the same scalar
C     product that SEAICE_FGMRES does
       CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)
       CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid)
       resLoc = SQRT(resLoc)
C     some output diagnostics
       IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(2A,2(1XI6),3E12.5)') 
     &       ' S/R SEAICE_JFNK_UPDATE: Newton iter, LSiter, ',
     &       'facLS, JFNKresidual, resLoc = ',
     &        newtonIter, l, facLS, JFNKresidual, resLoc
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &       SQUEEZE_RIGHT, myThid )
        _END_MASTER( myThid )
       ENDIF
C     Get ready for the next iteration: after adding du/vIce in the first
C     iteration, we substract 0.5*du/vIce from u/vIce in the next
C     iterations, 0.25*du/vIce in the second, etc.
       facLS = - 0.5 _d 0 * ABS(facLS)
      ENDDO
C     This is the new residual
      JFNKresidual = resLoc

#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */

      RETURN
      END