C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.31 2017/06/08 15:33:50 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
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
#ifdef SEAICE_ALLOW_JFNK
C !FUNCTIONS:
LOGICAL DIFFERENT_MULTIPLE
EXTERNAL
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 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
INTEGER iCode
_RL JFNKresidual
_RL JFNKresidualKm1
C parameters to compute convergence criterion
_RL JFNKgamma_lin
_RL FGMRESeps
_RL JFNKtol
C backward differences extrapolation factors
_RL bdfFac, bdfAlpha
C
_RL recip_deltaT
LOGICAL JFNKconverged, krylovConverged
LOGICAL writeNow
CHARACTER*(MAX_LEN_MBUF) msgBuf
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)
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
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 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
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
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
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*(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)*duIcNm1(I,J,bi,bj)*recip_deltaT
FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
& + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
ENDDO
ENDDO
CML ENDIF
ENDDO
ENDDO
C Start nonlinear Newton iteration: outer loop iteration
DO WHILE ( newtonIter.LT.SEAICEnonLinIterMax .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)
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)
ENDDO
ENDDO
ENDDO
ENDDO
C compute convergence criterion for linear preconditioned FGMRES
JFNKgamma_lin = JFNKgamma_lin_max
IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter
& .AND.JFNKresidual.LT.JFNKres_t ) THEN
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
C save the residual for the next iteration
JFNKresidualKm1 = JFNKresidual
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
JFNKconverged = JFNKresidual.LT.JFNKtol
& .OR.JFNKresidual.EQ.0. _d 0
C do Krylov loop only if convergence is not reached
IF ( .NOT.JFNKconverged ) THEN
C start Krylov iteration (FGMRES)
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 bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1,nVec
rhs(j,bi,bj) = - rhs(j,bi,bj)
ENDDO
ENDDO
ENDDO
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 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,SEAICElinearIterMax,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
IF (iCode.EQ.1) THEN
C Call preconditioner
IF ( SEAICEpreconLinIter .GT. 0 )
& CALL SEAICE_PRECONDITIONER(
U duIce, dvIce,
I zetaPre, etaPre, etaZpre, zetaZpre, 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 =
& SEAICEnonLinIterMax*(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.SEAICElinearIterMax ) 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=SEAICEnonLinTol*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-- Output diagnostics
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. SEAICEnonLinIterMax ) 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. SEAICEnonLinIterMax ) 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 resTmp :: vector version of the residuals
INTEGER nVec
PARAMETER ( nVec = 2*sNx*sNy )
_RL resTmp (nVec,1,nSx,nSy)
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 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 Determine, if we need more iterations
doLineSearch = resLoc .GE. JFNKresidual
C Limit the maximum number of iterations arbitrarily to four
doLineSearch = doLineSearch .AND. l .LT. 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 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_JFNK */
RETURN
END