pkg/seaice/seaice_jfnk.F

2fd913c523 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
a88fc20c6c Mart*0003 C--  File seaice_jfnk.F: seaice jfnk dynamical solver S/R:
                0004 C--   Contents
                0005 C--   o SEAICE_JFNK
                0006 C--   o SEAICE_JFNK_UPDATE
                0007
2fd913c523 Mart*0008 CBOP
                0009 C     !ROUTINE: SEAICE_JFNK
                0010 C     !INTERFACE:
                0011       SUBROUTINE SEAICE_JFNK( myTime, myIter, myThid )
                0012
                0013 C     !DESCRIPTION: \bv
                0014 C     *==========================================================*
a88fc20c6c Mart*0015 C     | SUBROUTINE SEAICE_JFNK
2fd913c523 Mart*0016 C     | o Ice dynamics using a Jacobian-free Newton-Krylov solver
                0017 C     |   following J.-F. Lemieux et al. Improving the numerical
                0018 C     |   convergence of viscous-plastic sea ice models with the
                0019 C     |   Jacobian-free Newton-Krylov method. J. Comp. Phys. 229,
                0020 C     |   2840-2852 (2010).
                0021 C     | o The logic follows JFs code.
                0022 C     *==========================================================*
                0023 C     | written by Martin Losch, Oct 2012
                0024 C     *==========================================================*
                0025 C     \ev
                0026
                0027 C     !USES:
                0028       IMPLICIT NONE
                0029
                0030 C     === Global variables ===
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
                0034 #include "DYNVARS.h"
                0035 #include "GRID.h"
                0036 #include "SEAICE_SIZE.h"
                0037 #include "SEAICE_PARAMS.h"
                0038 #include "SEAICE.h"
                0039
                0040 C     !INPUT/OUTPUT PARAMETERS:
                0041 C     === Routine arguments ===
                0042 C     myTime :: Simulation time
                0043 C     myIter :: Simulation timestep number
                0044 C     myThid :: my Thread Id. number
                0045       _RL     myTime
                0046       INTEGER myIter
                0047       INTEGER myThid
                0048
45315406aa Mart*0049 #if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_JFNK )
a711ed7dab Mart*0050 C     !FUNCTIONS:
                0051       LOGICAL  DIFFERENT_MULTIPLE
                0052       EXTERNAL DIFFERENT_MULTIPLE
2fd913c523 Mart*0053
740ca362e1 Mart*0054 C     !LOCAL VARIABLES:
                0055 C     === Local variables ===
79df32c3f1 Mart*0056 C     i,j,bi,bj :: loop indices
                0057       INTEGER i,j,bi,bj
2fd913c523 Mart*0058 C     loop indices
a711ed7dab Mart*0059       INTEGER newtonIter
                0060       INTEGER krylovIter, krylovFails
438b73e6d6 Mart*0061       INTEGER totalKrylovItersLoc, totalNewtonItersLoc
afd8ed70b7 Mart*0062 C     FGMRES parameters
                0063 C     im      :: size of Krylov space
                0064 C     ifgmres :: interation counter
                0065       INTEGER im
                0066       PARAMETER ( im = 50 )
                0067       INTEGER ifgmres
a711ed7dab Mart*0068 C     FGMRES flag that determines amount of output messages of fgmres
                0069       INTEGER iOutFGMRES
                0070 C     FGMRES flag that indicates what fgmres wants us to do next
2fd913c523 Mart*0071       INTEGER iCode
438b73e6d6 Mart*0072       _RL     JFNKresidual
2fd913c523 Mart*0073       _RL     JFNKresidualKm1
                0074 C     parameters to compute convergence criterion
1f3ad2d627 Mart*0075       _RL     JFNKgamma_lin
2fd913c523 Mart*0076       _RL     FGMRESeps
                0077       _RL     JFNKtol
e501eee760 Mart*0078 C     backward differences extrapolation factors
                0079       _RL bdfFac, bdfAlpha
6cbc659de0 Mart*0080 C
2fd913c523 Mart*0081       _RL     recip_deltaT
                0082       LOGICAL JFNKconverged, krylovConverged
f670f4bfe2 Mart*0083       LOGICAL writeNow
2fd913c523 Mart*0084       CHARACTER*(MAX_LEN_MBUF) msgBuf
10b2761ffa Jean*0085
2fd913c523 Mart*0086 C     u/vIceRes :: residual of sea-ice momentum equations
                0087       _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0088       _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e501eee760 Mart*0089 C     extra time level required for backward difference time stepping
6cbc659de0 Mart*0090       _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0091       _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2fd913c523 Mart*0092 C     du/vIce   :: ice velocity increment to be added to u/vIce
                0093       _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0094       _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
10b2761ffa Jean*0095 C     precomputed (= constant per Newton iteration) versions of
48c0fbd9c1 Mart*0096 C     zeta, eta, and DWATN, press
                0097       _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*0098       _RL zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48c0fbd9c1 Mart*0099       _RL etaPre  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f9529640cf Mart*0100       _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48c0fbd9c1 Mart*0101       _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
afd8ed70b7 Mart*0102 C     work arrays
                0103       _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy)
                0104       _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy)
                0105       _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy)
2fd913c523 Mart*0106 CEOP
                0107
                0108 C     Initialise
a711ed7dab Mart*0109       newtonIter          = 0
                0110       krylovFails         = 0
                0111       totalKrylovItersLoc = 0
                0112       JFNKconverged       = .FALSE.
                0113       JFNKtol             = 0. _d 0
                0114       JFNKresidual        = 0. _d 0
                0115       JFNKresidualKm1     = 0. _d 0
                0116       FGMRESeps           = 0. _d 0
                0117       recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn
                0118
                0119       iOutFGMRES=0
c09384f962 Mart*0120 C     with iOutFgmres=1, seaice_fgmres prints the residual at each iteration
                0121       IF ( debugLevel.GE.debLevC .AND.
a711ed7dab Mart*0122      &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
                0123      &     iOutFGMRES=1
                0124
e501eee760 Mart*0125 C     backward difference extrapolation factors
                0126       bdfFac = 0. _d 0
                0127       IF ( SEAICEuseBDF2 ) THEN
                0128        IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
                0129         bdfFac = 0. _d 0
6cbc659de0 Mart*0130        ELSE
e501eee760 Mart*0131         bdfFac = 0.5 _d 0
6cbc659de0 Mart*0132        ENDIF
                0133       ENDIF
c512e371cc drin*0134       bdfAlpha = 1. _d 0 + bdfFac
6cbc659de0 Mart*0135
2fd913c523 Mart*0136       DO bj=myByLo(myThid),myByHi(myThid)
                0137        DO bi=myBxLo(myThid),myBxHi(myThid)
10b2761ffa Jean*0138         DO J=1-OLy,sNy+OLy
                0139          DO I=1-OLx,sNx+OLx
2fd913c523 Mart*0140           uIceRes(I,J,bi,bj) = 0. _d 0
                0141           vIceRes(I,J,bi,bj) = 0. _d 0
                0142           duIce  (I,J,bi,bj) = 0. _d 0
                0143           dvIce  (I,J,bi,bj) = 0. _d 0
6cbc659de0 Mart*0144          ENDDO
                0145         ENDDO
                0146 C     cycle ice velocities
                0147         DO J=1-OLy,sNy+OLy
                0148          DO I=1-OLx,sNx+OLx
c512e371cc drin*0149           duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha
e501eee760 Mart*0150      &         + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac
c512e371cc drin*0151           dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha
e501eee760 Mart*0152      &         + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac
2fd913c523 Mart*0153           uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
                0154           vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
                0155          ENDDO
                0156         ENDDO
65e94703fa Mart*0157 C     As long as IMEX is not properly implemented leave this commented out
                0158 CML        IF ( .NOT.SEAICEuseIMEX ) THEN
2fd913c523 Mart*0159 C     Compute things that do no change during the Newton iteration:
10b2761ffa Jean*0160 C     sea-surface tilt and wind stress:
52a1d2474b Mart*0161 C     FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT
10b2761ffa Jean*0162         DO J=1-OLy,sNy+OLy
                0163          DO I=1-OLx,sNx+OLx
2fd913c523 Mart*0164           FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj)
6cbc659de0 Mart*0165      &         + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT
2fd913c523 Mart*0166           FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
6cbc659de0 Mart*0167      &         + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
2fd913c523 Mart*0168          ENDDO
                0169         ENDDO
65e94703fa Mart*0170 CML        ENDIF
2fd913c523 Mart*0171        ENDDO
                0172       ENDDO
                0173 C     Start nonlinear Newton iteration: outer loop iteration
79df32c3f1 Mart*0174       DO WHILE ( newtonIter.LT.SEAICEnonLinIterMax .AND.
2fd913c523 Mart*0175      &     .NOT.JFNKconverged )
                0176        newtonIter = newtonIter + 1
                0177 C     Compute initial residual F(u), (includes computation of global
                0178 C     variables DWATN, zeta, and eta)
10b2761ffa Jean*0179        IF ( newtonIter .EQ. 1 ) CALL SEAICE_JFNK_UPDATE(
                0180      I      duIce, dvIce,
a88fc20c6c Mart*0181      U      uIce, vIce, JFNKresidual,
                0182      O      uIceRes, vIceRes,
                0183      I      newtonIter, myTime, myIter, myThid )
2fd913c523 Mart*0184 C     local copies of precomputed coefficients that are to stay
                0185 C     constant for the preconditioner
                0186        DO bj=myByLo(myThid),myByHi(myThid)
                0187         DO bi=myBxLo(myThid),myBxHi(myThid)
10b2761ffa Jean*0188          DO j=1-OLy,sNy+OLy
                0189           DO i=1-OLx,sNx+OLx
bfe5e03139 Mart*0190            zetaPre(I,J,bi,bj) =  zeta(I,J,bi,bj)
0fd6b2de1a Mart*0191            zetaZPre(I,J,bi,bj)= zetaZ(I,J,bi,bj)
bfe5e03139 Mart*0192             etaPre(I,J,bi,bj) =   eta(I,J,bi,bj)
                0193            etaZPre(I,J,bi,bj) =  etaZ(I,J,bi,bj)
                0194            dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
2fd913c523 Mart*0195           ENDDO
                0196          ENDDO
                0197         ENDDO
                0198        ENDDO
                0199 C     compute convergence criterion for linear preconditioned FGMRES
                0200        JFNKgamma_lin = JFNKgamma_lin_max
f6f4a9e227 Mart*0201        IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter
2fd913c523 Mart*0202      &      .AND.JFNKresidual.LT.JFNKres_t ) THEN
1f3ad2d627 Mart*0203 C     Eisenstat and Walker (1996), eq.(2.6)
                0204         JFNKgamma_lin = SEAICE_JFNKphi
                0205      &       *( JFNKresidual/JFNKresidualKm1 )**SEAICE_JFNKalpha
2fd913c523 Mart*0206         JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin)
                0207         JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin)
                0208        ENDIF
                0209 C     save the residual for the next iteration
                0210        JFNKresidualKm1 = JFNKresidual
10b2761ffa Jean*0211
2fd913c523 Mart*0212 C     The Krylov iteration using FGMRES, the preconditioner is LSOR
                0213 C     for now. The code is adapted from SEAICE_LSR, but heavily stripped
                0214 C     down.
                0215 C     krylovIter is mapped into "its" in seaice_fgmres and is incremented
                0216 C     in that routine
                0217        krylovIter    = 0
                0218        iCode         = 0
10b2761ffa Jean*0219
2fd913c523 Mart*0220        JFNKconverged = JFNKresidual.LT.JFNKtol
0d02be7d13 Mart*0221      &      .OR.JFNKresidual.EQ.0. _d 0
10b2761ffa Jean*0222
2fd913c523 Mart*0223 C     do Krylov loop only if convergence is not reached
10b2761ffa Jean*0224
2fd913c523 Mart*0225        IF ( .NOT.JFNKconverged ) THEN
10b2761ffa Jean*0226
2fd913c523 Mart*0227 C     start Krylov iteration (FGMRES)
10b2761ffa Jean*0228
2fd913c523 Mart*0229         krylovConverged = .FALSE.
                0230         FGMRESeps = JFNKgamma_lin * JFNKresidual
afd8ed70b7 Mart*0231 C     map first guess sol; it is zero because the solution is a correction
                0232        CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid)
                0233 C     map rhs and change its sign because we are solving J*u = -F
79df32c3f1 Mart*0234         CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,rhs,.TRUE.,myThid)
                0235         DO bj=myByLo(myThid),myByHi(myThid)
                0236          DO bi=myBxLo(myThid),myBxHi(myThid)
                0237           DO j=1,nVec
                0238            rhs(j,bi,bj) = - rhs(j,bi,bj)
                0239           ENDDO
                0240          ENDDO
                0241         ENDDO
10b2761ffa Jean*0242         DO WHILE ( .NOT.krylovConverged )
2fd913c523 Mart*0243 C     solution vector sol = du/vIce
                0244 C     residual vector (rhs) Fu = u/vIceRes
10b2761ffa Jean*0245 C     output work vectors wk1, -> input work vector wk2
                0246
afd8ed70b7 Mart*0247 C     map preconditioner results or Jacobian times vector,
                0248 C     stored in du/vIce to wk2, for iCode=0, wk2 is set to zero,
                0249 C     because du/vIce = 0
                0250          CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid)
c512e371cc drin*0251
afd8ed70b7 Mart*0252          CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,
                0253      U        vv,w,wk1,wk2,
79df32c3f1 Mart*0254      I        FGMRESeps,SEAICElinearIterMax,iOutFGMRES,
afd8ed70b7 Mart*0255      U        iCode,
                0256      I        myThid)
c512e371cc drin*0257
afd8ed70b7 Mart*0258          IF ( iCode .EQ. 0 ) THEN
                0259 C     map sol(ution) vector to du/vIce
                0260           CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid)
                0261          ELSE
                0262 C     map work vector to du/vIce to either compute a preconditioner
                0263 C     solution (wk1=rhs) or a Jacobian times wk1
                0264           CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid)
                0265          ENDIF
                0266 C     Fill overlaps in updated fields
                0267          CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)
2fd913c523 Mart*0268 C     FGMRES returns iCode either asking for an new preconditioned vector
                0269 C     or product of matrix (Jacobian) times vector. For iCode = 0, terminate
                0270 C     iteration
                0271          IF (iCode.EQ.1) THEN
10b2761ffa Jean*0272 C     Call preconditioner
79df32c3f1 Mart*0273           IF ( SEAICEpreconLinIter .GT. 0 )
10b2761ffa Jean*0274      &         CALL SEAICE_PRECONDITIONER(
                0275      U         duIce, dvIce,
0fd6b2de1a Mart*0276      I         zetaPre, etaPre, etaZpre, zetaZpre, dwatPre,
2fd913c523 Mart*0277      I         newtonIter, krylovIter, myTime, myIter, myThid )
                0278          ELSEIF (iCode.GE.2) THEN
                0279 C     Compute Jacobian times vector
                0280           CALL SEAICE_JACVEC(
                0281      I         uIce, vIce, uIceRes, vIceRes,
10b2761ffa Jean*0282      U         duIce, dvIce,
2fd913c523 Mart*0283      I         newtonIter, krylovIter, myTime, myIter, myThid )
                0284          ENDIF
                0285          krylovConverged = iCode.EQ.0
                0286 C     End of Krylov iterate
                0287         ENDDO
a711ed7dab Mart*0288         totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
2fd913c523 Mart*0289 C     some output diagnostics
                0290         IF ( debugLevel.GE.debLevA ) THEN
a711ed7dab Mart*0291          _BEGIN_MASTER( myThid )
10b2761ffa Jean*0292          totalNewtonItersLoc =
79df32c3f1 Mart*0293      &        SEAICEnonLinIterMax*(myIter-nIter0)+newtonIter
10b2761ffa Jean*0294          WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
438b73e6d6 Mart*0295      &        ' S/R SEAICE_JFNK: Newton iterate / total, ',
                0296      &        'JFNKgamma_lin, initial norm = ',
                0297      &        newtonIter, totalNewtonItersLoc,
                0298      &        JFNKgamma_lin,JFNKresidual
                0299          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0300      &        SQUEEZE_RIGHT, myThid )
2fd913c523 Mart*0301          WRITE(msgBuf,'(3(A,I6))')
10b2761ffa Jean*0302      &        ' S/R SEAICE_JFNK: Newton iterate / total = ',newtonIter,
438b73e6d6 Mart*0303      &        ' / ', totalNewtonItersLoc,
2fd913c523 Mart*0304      &        ', Nb. of FGMRES iterations = ', krylovIter
                0305          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0306      &        SQUEEZE_RIGHT, myThid )
a711ed7dab Mart*0307          _END_MASTER( myThid )
2fd913c523 Mart*0308         ENDIF
79df32c3f1 Mart*0309         IF ( krylovIter.EQ.SEAICElinearIterMax ) THEN
a711ed7dab Mart*0310          krylovFails = krylovFails + 1
2fd913c523 Mart*0311         ENDIF
2e75855dde Mart*0312 C     Set the stopping criterion for the Newton iteration and the
                0313 C     criterion for the transition from accurate to approximate FGMRES
                0314         IF ( newtonIter .EQ. 1 ) THEN
79df32c3f1 Mart*0315          JFNKtol=SEAICEnonLinTol*JFNKresidual
2e75855dde Mart*0316          IF ( JFNKres_tFac .NE. UNSET_RL )
                0317      &        JFNKres_t = JFNKresidual * JFNKres_tFac
                0318         ENDIF
2fd913c523 Mart*0319 C     Update linear solution vector and return to Newton iteration
a88fc20c6c Mart*0320 C     Do a linesearch if necessary, and compute a new residual.
                0321 C     Note that it should be possible to do the following operations
                0322 C     at the beginning of the Newton iteration, thereby saving us from
                0323 C     the extra call of seaice_jfnk_update, but unfortunately that
                0324 C     changes the results, so we leave the stuff here for now.
10b2761ffa Jean*0325         CALL SEAICE_JFNK_UPDATE(
                0326      I       duIce, dvIce,
a88fc20c6c Mart*0327      U       uIce, vIce, JFNKresidual,
                0328      O       uIceRes, vIceRes,
                0329      I       newtonIter, myTime, myIter, myThid )
                0330 C     reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
2fd913c523 Mart*0331         DO bj=myByLo(myThid),myByHi(myThid)
                0332          DO bi=myBxLo(myThid),myBxHi(myThid)
10b2761ffa Jean*0333           DO J=1-OLy,sNy+OLy
                0334            DO I=1-OLx,sNx+OLx
68c9371d34 Mart*0335             duIce(I,J,bi,bj)= 0. _d 0
                0336             dvIce(I,J,bi,bj)= 0. _d 0
2fd913c523 Mart*0337            ENDDO
                0338           ENDDO
                0339          ENDDO
                0340         ENDDO
                0341        ENDIF
                0342 C     end of Newton iterate
                0343       ENDDO
10b2761ffa Jean*0344
a711ed7dab Mart*0345 C--   Output diagnostics
10b2761ffa Jean*0346
186a000033 Mart*0347       IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
c512e371cc drin*0348 C     Only Master Thread updates counters in common block:
                0349       _BEGIN_MASTER(myThid)
                0350 C     Cumulate some diagnostic counters for the JFNK solver
186a000033 Mart*0351        totalJFNKtimeSteps = totalJFNKtimeSteps + 1
                0352        totalNewtonIters   = totalNewtonIters + newtonIter
                0353        totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
2fd913c523 Mart*0354 C     Record failure
186a000033 Mart*0355        totalKrylovFails   = totalKrylovFails + krylovFails
79df32c3f1 Mart*0356        IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN
10b2761ffa Jean*0357         totalNewtonFails = totalNewtonFails + 1
186a000033 Mart*0358        ENDIF
c512e371cc drin*0359        _END_MASTER( myThid )
a711ed7dab Mart*0360       ENDIF
                0361 C     Decide whether it is time to dump and reset the counter
f670f4bfe2 Mart*0362       writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
10b2761ffa Jean*0363      &     myTime+deltaTClock, deltaTClock)
f670f4bfe2 Mart*0364 #ifdef ALLOW_CAL
                0365       IF ( useCAL ) THEN
10b2761ffa Jean*0366        CALL CAL_TIME2DUMP(
f670f4bfe2 Mart*0367      I      zeroRL, SEAICE_monFreq,  deltaTClock,
                0368      U      writeNow,
c512e371cc drin*0369      I      myTime+deltaTClock, myIter+1, myThid )
f670f4bfe2 Mart*0370       ENDIF
                0371 #endif
                0372       IF ( writeNow ) THEN
a711ed7dab Mart*0373        _BEGIN_MASTER( myThid )
10b2761ffa Jean*0374        WRITE(msgBuf,'(A)')
a711ed7dab Mart*0375      &' // ======================================================='
                0376        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0377      &      SQUEEZE_RIGHT, myThid )
                0378        WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
                0379        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0380      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0381        WRITE(msgBuf,'(A)')
a711ed7dab Mart*0382      &' // ======================================================='
                0383        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0384      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0385        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0386      &      ' %JFNK_MON: time step              = ', myIter+1
                0387        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0388      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0389        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0390      &      ' %JFNK_MON: Nb. of time steps      = ', totalJFNKtimeSteps
                0391        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0392      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0393        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0394      &      ' %JFNK_MON: Nb. of Newton steps    = ', totalNewtonIters
                0395        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0396      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0397        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0398      &      ' %JFNK_MON: Nb. of Krylov steps    = ', totalKrylovIters
                0399        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0400      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0401        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0402      &      ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
                0403        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0404      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0405        WRITE(msgBuf,'(A,I10)')
a711ed7dab Mart*0406      &      ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
                0407        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0408      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0409        WRITE(msgBuf,'(A)')
a711ed7dab Mart*0410      &' // ======================================================='
                0411        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0412      &      SQUEEZE_RIGHT, myThid )
55c264c0e4 Mart*0413        WRITE(msgBuf,'(A)') ' // End JFNK statistics'
a711ed7dab Mart*0414        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0415      &      SQUEEZE_RIGHT, myThid )
10b2761ffa Jean*0416        WRITE(msgBuf,'(A)')
a711ed7dab Mart*0417      &' // ======================================================='
                0418        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0419      &      SQUEEZE_RIGHT, myThid )
c512e371cc drin*0420 C     Reset and start again
a711ed7dab Mart*0421        totalJFNKtimeSteps = 0
                0422        totalNewtonIters   = 0
                0423        totalKrylovIters   = 0
                0424        totalKrylovFails   = 0
                0425        totalNewtonFails   = 0
c512e371cc drin*0426        _END_MASTER( myThid )
a711ed7dab Mart*0427       ENDIF
                0428
                0429 C     Print more debugging information
                0430       IF ( debugLevel.GE.debLevA ) THEN
79df32c3f1 Mart*0431        IF ( newtonIter .EQ. SEAICEnonLinIterMax ) THEN
a711ed7dab Mart*0432         _BEGIN_MASTER( myThid )
10b2761ffa Jean*0433         WRITE(msgBuf,'(A,I10)')
2fd913c523 Mart*0434      &       ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
a711ed7dab Mart*0435      &       myIter+1
2fd913c523 Mart*0436         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0437      &       SQUEEZE_RIGHT, myThid )
a711ed7dab Mart*0438         _END_MASTER( myThid )
2fd913c523 Mart*0439        ENDIF
a711ed7dab Mart*0440        IF ( krylovFails .GT. 0 ) THEN
                0441         _BEGIN_MASTER( myThid )
10b2761ffa Jean*0442         WRITE(msgBuf,'(A,I4,A,I10)')
2fd913c523 Mart*0443      &       ' S/R SEAICE_JFNK: FGMRES did not converge ',
a711ed7dab Mart*0444      &       krylovFails, ' times in timestep ', myIter+1
2fd913c523 Mart*0445         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0446      &       SQUEEZE_RIGHT, myThid )
a711ed7dab Mart*0447         _END_MASTER( myThid )
2fd913c523 Mart*0448        ENDIF
a711ed7dab Mart*0449        _BEGIN_MASTER( myThid )
10b2761ffa Jean*0450        WRITE(msgBuf,'(A,I6,A,I10)')
2fd913c523 Mart*0451      &      ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
a711ed7dab Mart*0452      &      totalKrylovItersLoc, ' in timestep ', myIter+1
                0453        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0454      &      SQUEEZE_RIGHT, myThid )
                0455        _END_MASTER( myThid )
2fd913c523 Mart*0456       ENDIF
                0457
a88fc20c6c Mart*0458       RETURN
                0459       END
                0460
740ca362e1 Mart*0461 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
a88fc20c6c Mart*0462 CBOP
                0463 C     !ROUTINE: SEAICE_JFNK_UPDATE
                0464 C     !INTERFACE:
                0465
10b2761ffa Jean*0466       SUBROUTINE SEAICE_JFNK_UPDATE(
                0467      I     duIce, dvIce,
a88fc20c6c Mart*0468      U     uIce, vIce, JFNKresidual,
                0469      O     uIceRes, vIceRes,
                0470      I     newtonIter, myTime, myIter, myThid )
                0471
                0472 C     !DESCRIPTION: \bv
                0473 C     *==========================================================*
                0474 C     | SUBROUTINE SEAICE_JFNK_UPDATE
                0475 C     | o Update velocities with incremental solutions of FGMRES
                0476 C     | o compute residual of updated solutions and do
                0477 C     | o linesearch:
                0478 C     |   reduce update until residual is smaller than previous
                0479 C     |   one (input)
                0480 C     *==========================================================*
                0481 C     | written by Martin Losch, Jan 2013
                0482 C     *==========================================================*
                0483 C     \ev
                0484
                0485 C     !USES:
                0486       IMPLICIT NONE
                0487
                0488 C     === Global variables ===
                0489 #include "SIZE.h"
                0490 #include "EEPARAMS.h"
                0491 #include "PARAMS.h"
                0492 #include "SEAICE_SIZE.h"
                0493 #include "SEAICE_PARAMS.h"
                0494
                0495 C     !INPUT/OUTPUT PARAMETERS:
                0496 C     === Routine arguments ===
                0497 C     myTime :: Simulation time
                0498 C     myIter :: Simulation timestep number
                0499 C     myThid :: my Thread Id. number
                0500 C     newtonIter :: current iterate of Newton iteration
                0501       _RL     myTime
                0502       INTEGER myIter
                0503       INTEGER myThid
                0504       INTEGER newtonIter
                0505 C     JFNKresidual :: Residual at the beginning of the FGMRES iteration,
                0506 C                     changes with newtonIter (updated)
                0507       _RL     JFNKresidual
                0508 C     du/vIce   :: ice velocity increment to be added to u/vIce (input)
                0509       _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0510       _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0511 C     u/vIce    :: ice velocity increment to be added to u/vIce (updated)
                0512       _RL uIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0513       _RL vIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0514 C     u/vIceRes :: residual of sea-ice momentum equations (output)
                0515       _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0516       _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0517
740ca362e1 Mart*0518 C     !LOCAL VARIABLES:
                0519 C     === Local variables ===
a88fc20c6c Mart*0520 C     i,j,bi,bj :: loop indices
                0521       INTEGER i,j,bi,bj
                0522       INTEGER l
                0523       _RL     resLoc, facLS
                0524       LOGICAL doLineSearch
                0525 C     nVec    :: size of the input vector(s)
10b2761ffa Jean*0526 C     resTmp  :: vector version of the residuals
a88fc20c6c Mart*0527       INTEGER nVec
                0528       PARAMETER ( nVec  = 2*sNx*sNy )
                0529       _RL resTmp (nVec,1,nSx,nSy)
10b2761ffa Jean*0530
a88fc20c6c Mart*0531       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0532 CEOP
                0533
                0534 C     Initialise some local variables
                0535       l = 0
                0536       resLoc = JFNKresidual
                0537       facLS = 1. _d 0
                0538       doLineSearch = .TRUE.
                0539       DO WHILE ( doLineSearch )
c704c5a1ef Mart*0540 C     The computation of the line search factor is complicated because
                0541 C     we overwrite the solution, so that we have to reconstruct what
                0542 C     should be just u_l = u_0 + (1-gamma)**l *du. Instead, after adding
                0543 C     du/vIce in the first iteration (l=0), we substract a decreasing
                0544 C     part of du/vIce from u/vIce in the next iterations until the
                0545 C     residual is smaller than the initial residual.
                0546        IF ( l .GT. 0 ) facLS = - SEAICE_JFNK_lsGamma
                0547      &       *(1. _d 0 - SEAICE_JFNK_lsGamma)**(l-1)
a88fc20c6c Mart*0548 C     Create update
                0549        DO bj=myByLo(myThid),myByHi(myThid)
                0550         DO bi=myBxLo(myThid),myBxHi(myThid)
10b2761ffa Jean*0551          DO J=1-OLy,sNy+OLy
                0552           DO I=1-OLx,sNx+OLx
a88fc20c6c Mart*0553            uIce(I,J,bi,bj) = uIce(I,J,bi,bj)+facLS*duIce(I,J,bi,bj)
                0554            vIce(I,J,bi,bj) = vIce(I,J,bi,bj)+facLS*dvIce(I,J,bi,bj)
                0555           ENDDO
                0556          ENDDO
                0557         ENDDO
                0558        ENDDO
                0559 C     Compute current residual F(u), (includes re-computation of global
                0560 C     variables DWATN, zeta, and eta, i.e. they are different after this)
10b2761ffa Jean*0561        CALL SEAICE_CALC_RESIDUAL(
                0562      I      uIce, vIce,
                0563      O      uIceRes, vIceRes,
a88fc20c6c Mart*0564      I      newtonIter, 0, myTime, myIter, myThid )
                0565 C     Important: Compute the norm of the residual using the same scalar
                0566 C     product that SEAICE_FGMRES does
                0567        CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)
                0568        CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid)
                0569        resLoc = SQRT(resLoc)
175055f98c Mart*0570 C     Determine, if we need more iterations
10b2761ffa Jean*0571        doLineSearch = resLoc .GE. JFNKresidual
c704c5a1ef Mart*0572 C     Limit the maximum number of iterations arbitrarily to
                0573 C     SEAICE_JFNK_lsLmax (default is 4)
                0574        doLineSearch = doLineSearch .AND. l .LE. SEAICE_JFNK_lsLmax
175055f98c Mart*0575 C     For the first iteration du/vIce = 0 and there will be no
                0576 C     improvement of the residual possible, so we do only the first
                0577 C     iteration
                0578        IF ( newtonIter .EQ. 1 ) doLineSearch = .FALSE.
                0579 C     Only start a linesearch after some Newton iterations
                0580        IF ( newtonIter .LE. SEAICE_JFNK_lsIter ) doLineSearch = .FALSE.
a88fc20c6c Mart*0581 C     some output diagnostics
                0582        IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN
                0583         _BEGIN_MASTER( myThid )
c704c5a1ef Mart*0584         WRITE(msgBuf,'(2A,2(1XI6),1X,3E12.5)')
a88fc20c6c Mart*0585      &       ' S/R SEAICE_JFNK_UPDATE: Newton iter, LSiter, ',
                0586      &       'facLS, JFNKresidual, resLoc = ',
                0587      &        newtonIter, l, facLS, JFNKresidual, resLoc
                0588         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0589      &       SQUEEZE_RIGHT, myThid )
                0590         _END_MASTER( myThid )
                0591        ENDIF
c704c5a1ef Mart*0592 C     Increment counter
                0593        l = l + 1
a88fc20c6c Mart*0594       ENDDO
                0595 C     This is the new residual
                0596       JFNKresidual = resLoc
                0597
45315406aa Mart*0598 #endif /* SEAICE_CGRID and SEAICE_ALLOW_JFNK */
2fd913c523 Mart*0599
                0600       RETURN
                0601       END